Mechanical response of packings of non-spherical particles: A case study of 2D packings of circulo-lines

Granular materials represent fascinating examples of nonequilibrium physical systems that display complex, collective behavior [1, 2], such as stick-slip motion [3], shear banding [4], and segregation [5]. They are composed of discrete, macroscopic grains that interact via dissipative, frictional contact interactions. Granular materials occur in many important contexts, including numerous geological processes [6, 7], food and consumer product processing [8, 9], and robotics applications [10, 11]. Because granular media are highly dissipative, the grains do not move unless they are driven by gravity, fluid flow, or the system boundaries. When granular systems are compressed to sufficiently large packing fractions, they become jammed and possess a nonzero static shear modulus and other soid-like properties [12].

While most dry granular materials are composed of frictional, non-spherical grains, numerous computational and theoretical studies of dry granular packings have focused on the soft-particle model in which frictionless, spherical particles interact via the pairwise, purely repulsive potential energy [13]: $U(r_{ij})\propto(1-r_{ij}/\sigma_{ij})^{\alpha}\Theta(1-r_{ij}/\sigma_{ij})$, where $r_{ij}$ is the separation between the centers of mass of particles $i$ and $j$, $\sigma_{ij}$ is their average diameter, the exponents $\alpha=2$ and $5/2$ are set for purely repulsive linear and Hertzian spring interactions [14], and $\Theta(.)$ is the Heaviside step function that ensures the potential energy is nonzero only when the particles are in contact. For frictionless, spherical particles, the jamming transition occurs when the number of interparticle contacts $N_{c}$ equals or exceeds the isostatic value $N_{c}^{0}=dN-d+1$, where $d$ is the spatial dimension and $N$ is the number of non-rattler particles [15]. When the system is compressed above jamming onset, the pressure increases from zero and the jammed packing develops nonzero bulk and shear moduli.

A hallmark of the jamming transition in static packings of frictionless, spherical particles is that the ensemble-averaged contact number and shear modulus $\langle G\rangle$ scale as a power-law in the pressure $P$ [12, 16], above a characteristic pressure $P^{**}$ that decreases with increasing system size, as the packings are isotropically compressed above jamming onset [17]. For example, $\langle G\rangle\sim P^{1/2}$ for $P>P^{**}$ for packings of spherical particles with purely repulsive, linear spring interactions. Several previous experimental studies of compressed emulsions [18, 19] and packings of thin granular cylinders [20] have also found power-law scaling of the contact number and shear modulus with pressure during isotropic compression. In previous computational studies of packings of frictionless, spherical particles [21], we showed that there are two important contributions to the ensemble-averaged shear modulus: $G_{f}$ from geometrical families and $G_{r}$ from changes in the interparticle contact network during compression. For isotropic compression, jammed packings within a geometrical family are mechanically stable packings with different pressures that are related to each other by continuous, quasistatic changes in packing fraction with no changes in the interparticle contact network. $G_{r}$ includes discontinuities in the shear modulus that arise from point and jump changes in the contact network [22, 23]. Point changes involve the addition or removal of a single interparticle contact (or multiple contacts when a rattler particle is added or removed from the contact network) without significant particle motion. Jump changes are caused by mechanical instabilities and typically involve multiple changes in the contact network and collective particle motion.

For frictionless, spherical particles, the shear modulus along a geometrical family decreases with increasing pressure during isotropic compression (and no shearing of the boundaries). For systems with purely repulsive linear spring interactions, $G_{f}$ decreases roughly linearly with pressure, $G_{f}(P)/G_{0}\sim 1-P/P_{0}$, where $G_{0}$ is the shear modulus at $P=0$ and $P_{0}$ is the pressure at which $G_{f}(P)=0$, although there are deviations from this simple form as the system approaches contact network changes [24, 21]. Thus, the increase in the ensemble-averaged shear modulus, $\langle G\rangle$, with pressure for packings of spherical particles arises from $G_{r}$, where changes in the contact network cause discontinuous upward jumps in the shear modulus. We will show below that the sign of the second, pressure-dependent term in $G_{f}$ (i.e. whether $G_{f}$ increases or decreases with pressure) is determined by the curvature of geometrical families in the strain direction.

Several prior computational studies have shown that the form of $\langle G(P)\rangle$ for jammed packings depends on the particle shape [25, 26]. For example, the pressure-depenent shear modulus for packings of frictionless ellipse-shaped particles with repulsive linear spring interactions scales as $\langle G(P)\rangle\sim P^{\beta}$, where $\beta\sim 1$ over a wide range of pressure, $P^{**}<P<P^{*}$, $P^{**}\sim N^{-2}$, and $P^{*}$ does not depend strongly on system size. For $P>P^{*}$, the power-law scaling crosses over to $\langle G\rangle\sim P^{1/2}$, as found for jammed packings of spherical particles [17]. However, the ensemble-averaged shear modulus for jammed packings of dimer-shaped particles (and other composite particles formed from bonded spherical particles) scales as $\langle G\rangle\sim P^{1/2}$ over the same range of pressures and for the same aspect ratios as those studied for packings of smooth, ellipse-shaped particles [26]. Since packings of frictionless ellipse-shaped particles are hypostatic at jamming onset, while packings of dimer-shaped particles are isostatic at jamming onset, it is possible that the change in the power-law scaling exponent from $\beta=0.5$ to $\sim 1$ is related to the presence of low-frequency quartic modes in the vibrational response for packings of ellipse-shaped and other non-spherical particles [27].

In this work, we address the question of what determines the power-law scaling exponent $\beta$ for packings of non-spherical particles. Does $G_{f}$ decrease with pressure for packings of non-spherical particles? Are the frequency of contact network changes and the magnitude of the discontinuities in the shear modulus at contact changes different from those for packings of spherical particles? As a case study, we investigate the pressure-dependent mechanical response of jammed packings of circulo-lines interacting via purely repulsive, linear spring interactions in two spatial dimensions (2D) as a function of aspect ratio.

We find several key results for the mechanical response of jammed packings of circulo-lines. First, the curvature of the variation of packing fraction with strain for geometrical families can be either negative or positive, and thus the shear modulus of geometrical families can either increase or decrease with pressure: $G_{f}/G_{0}\sim 1\pm P/P_{0}$ to linear order in pressure. We derive an exact expression for the pressure-dependent shear modulus of jammed packings and show that near jamming onset it can be approximated as $G_{f}\sim-\frac{1}{\phi}\left(\frac{dP}{d\gamma}\right)_{\phi}\left(\frac{d\phi}{d\gamma}\right)_{P}-\frac{P}{\phi}\left(\frac{\left(\frac{d\phi}{d\gamma}\right)_{P}}{d\gamma}\right)_{\phi}$ [28], where $\gamma$ is the shear strain. The first term tends to a constant, $G_{0}>0$, in the zero-pressure limit and the sign of the coefficient of the second term (that is roughly linear in $P$) is determined by the curvature of geometric families in the $\phi$-$\gamma$ plane. Second, we decompose $G_{f}=G_{a}-G_{na}$ for each first, low-pressure geometrical family into its affine and non-affine contributions. $G_{a}$ gives the response of the system to a globally affine change of the particle positions and system boundaries, while $G_{na}$ also includes particle motion in response to potential energy minimization. We find that the non-affine term plays an important role in determining $G_{f}$. In particular, the non-affine contribution can cause $G_{f}$ to increase with pressure, which does not occur in jammed packings of spherical particles. We also calculate the ensemble-averaged shear modulus $\langle G\rangle$ versus $P$ for jammed packings of circulo-lines over a range of aspect ratios ${\cal R}$ and system sizes. For packings of circulo-lines, we find that $\langle G\rangle\sim P^{\beta}$, where $\beta\sim 0.8$-$0.9$, over a range of pressures $P^{**}<P<P^{*}$, where $P^{**}\sim N^{-2}$ and $P^{*}$ decreases as ${\cal R}\rightarrow 1$ and does not depend strongly on system size. We find that the finite fraction of geometrical families with negative curvature in the $\phi$-$\gamma$ plane causes the power-law exponent $\beta$ for $\langle G\rangle$ versus $P$ to increase compared to that for jammed packings of spherical particles.

The remainder of the article is organized as follows. In Sec. II, we describe the interparticle potential and types of contacts that occur between pairs of circulo-lines, the numerical methods used to generate jammed packings of circulo-lines, and formulas for calculating the pressure, shear stress, and shear modulus. In Sec. III, we describe the results concerning geometrical families, non-affine contributions to the shear modulus, point and jump changes, and the ensemble-averaged shear modulus for jammed packings of circulo-lines. Finally, in Sec. IV, we provide the conclusions and point to promising future research directions. We also include an appendix that gives explicit expressions for the affine shear modulus for packings of circulo-lines with purely repulsive, linear spring interactions.

A circulo-line is the set of points that are equally distant from a line segment; it is thus a 2D shape composed of a rectangular middle region capped by two semi-circles on both ends. Fig. 1 (a) shows a circulo-line (labeled $i$) with diameter of the semi-circles, $\sigma_{i}$, and length of the middle line segment, $2l_{i}$. We will refer to the end points of the middle line segment as the foci. The aspect ratio of a circulo-line is $\mathcal{R}=(\sigma_{i}+2l_{i})/\sigma_{i}$, and ${\cal R}=1$ in the limit that the particle becomes coincident semi-circles. The asphericity of a circulo-line is $\mathcal{A}=p_{i}^{2}/4\pi a_{i}=\frac{(2\pi+4(\mathcal{R}_{i}-1))^{2}}{4\pi(\pi+4(\mathcal{R}_{i}-1))}$, where $p_{i}$ is the perimeter and $a_{i}$ is the area of the circulo-line. In these studies, we vary ${\cal R}-1$ and ${\cal A}-1$ over the ranges $10^{-2}$ to $2$ and $10^{-5}$ to $0.5$, respectively. We study bidisperse packings of circulo-lines to inhibit positional and orientational order. We consider packings in rhombic boxes with edge length $L$, periodic boundary conditions in the $x$- and $y$-directions, and $N/2$ large and $N/2$ small particles with end cap diameter ratio $\sigma_{l}/\sigma_{s}=1.4$, but the same mass $m$ and aspect ratio ${\cal R}$. To investigate the effect of system size, we studied packings with $N=64$, $128$, $256$, and $512$.

In Fig. 1 (b)-(d), we show that there are three ways in which two circulo-lines can make contact (or make small overlaps) with each other: end-end, end-middle, and middle-middle contacts. Fig. 1 (b) shows an end-end contact between circulo-lines $i$ and $j$. In this case, the relevant separation $r_{ij}$ between the circulo-lines is the separation between the two closest foci on $i$ and $j$. Another important distance is the separation $\lambda_{i}$ between the center of circulo-line $i$ and the point at which the line from the closest foci on $j$ that is perpendicular to the axis of circulo-line $i$ intersects the axis of $i$. For an end-end contact, the following three conditions must be satisfied: $r_{ij}<\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$, $\lambda_{i}>l_{i}$, and $\lambda_{j}>l_{j}$.

Fig. 1 (c) shows an end-middle contact where the end of circulo-line $j$ is in contact with the middle region of circulo-line $i$. For an end-middle contact, the relevant separation is the distance $r_{ji}$ between the closest focus on $i$ and the axis of $j$. Similarly, we can define the separation $r_{ij}$, but $r_{ji}\neq r_{ij}$. For an end-middle contact, the following four conditions must be satisfied: $r_{ij}<\sigma_{ij}$, $r_{ji}>\sigma_{ij}$, $\lambda_{i}<l_{i}$, and $\lambda_{j}>l_{j}$.

Fig. 1 (d) shows two circulo-lines for which the two middle regions are in contact. We model middle-middle contacts as two end-middle contacts, so that when a middle-middle contact becomes an end-end contact (from Fig. 1 (d) to (b)) or end-middle contact (from Fig. 1 (d) to (c)) due to small changes in the positions and orientations of the circulo-lines, the interparticle potential energy and forces are continuous. Thus, for a middle-middle contact, the following four conditions must be satisfied: $r_{ij}<\sigma_{ij}$, $r_{ji}<\sigma_{ij}$, $\lambda_{i}<l_{i}$, and $\lambda_{j}<l_{j}$.

For all three types of contacts, we use the pairwise, purely repulsive, linear spring potential

where $\epsilon$ is the characteristic energy scale of the interaction and $\Theta(.)$ is the Heaviside step function, ensuring that pairs of circulo-lines do not interact when then are not in contact. The total potential energy, $U=\sum_{i>j}U(r_{ij})$, and pair force, ${\vec{f}}_{ij}=(dU/dr_{ij}){\hat{r}}_{ij}$, are continuous as a function of the coordinates ${\vec{r}}_{i}$ of the centers of mass and orientations $\theta_{i}$ of the circulo-lines. We use $\epsilon$, $\epsilon/\sigma_{s}$, and $\epsilon/\sigma^{2}_{s}$ for the units of energy, force, and stress.

We employ the Love expression [29] to calculate the stress tensor:

where $f_{ij\alpha}$ is the $\alpha$-component of the force on particle $i$ due to particle $j$ and $\rho_{ij\beta}$ is the $\beta$-component of the vector pointing from the effective contact point between two overlapping circulo-lines $i$ and $j$ to the center of $j$. The effective contact point $C_{ij}$ is the midpoint of $r_{ij}$ as shown in Fig. 1 (e). We define the pressure as $P=(\Sigma_{xx}+\Sigma_{yy})/2$ and the shear stress as $\Sigma=-\Sigma_{xy}$.

To calculate the shear modulus, we first apply an affine simple shear strain step $\delta\gamma$, which changes the circulo-line center of mass positions and orientations,

and

together with Lees-Edwards periodic boundary conditions. ($x_{i}$,$y_{i}$) gives the position of the center of mass of the $i$th circulo-line and $\theta_{i}$ gives the angle that the axis of the circulo-line makes with the $x$-axis. After each simple shear strain step, we minimize the total potential energy using the FIRE algorithm and measure the shear stress $\Sigma$. We can then determine the shear modulus by calculating $G=d\Sigma/d\gamma$.

We employ an athermal, quasistatic isotropic compression protocol to generate packings of circulo-lines at jamming onset. We initialize the system with random positions and orientations of the circulo-lines in the dilute limit with packing fraction $\phi<10^{-2}$. We then compress the system by $\Delta\phi/\phi=2\times 10^{-3}$, minimize the total potential energy using the FIRE algorithm, and measure the pressure $P$. If $P<P_{t}$, where $P_{t}$ is the target pressure, we again compress the system by $\Delta\phi$ and minimize the total potential energy. If $P>P_{t}$, we return to the previous configuration and decrease the packing fraction increment by a factor of two. We continue this compression process until $|P-P_{t}|/P_{t}<10^{-5}$. For $P_{t}<10^{-7}$, we find that the number of interparticle contacts satisfies the effective isostatic condition: $N_{c}=N_{c}^{0}-N_{q}$, where $N_{c}^{0}=d_{f}N^{\prime}-d+1$, $d=2$ is the spatial dimension, the number of degrees of freedom per particle, $d_{f}=3$, $N^{\prime}=N-N_{r}-N_{s}/3$, $N_{r}$ is the number of rattler particles with too few contacts to constrain the $d_{f}$ degrees of freedom, $N_{s}$ is the number of slider particles with one unconstrained translational degree of freedom, and $N_{q}$ is the number of quartic modes of the dynamical matrix [27]. To investigate the pressure-dependent mechanical properties, we start with systems at $P_{t}=10^{-7}$ and then successively increase the target pressure over the range $10^{-7}<P_{t}<10^{-2}$.

Here, we describe our main results in five subsections. In Sec. III.1, we generalize the concept of geometrical families of jammed packings in the $\phi$-$\gamma$ plane to packings of circulo-lines and show that the curvature $d^{2}\phi/d\gamma^{2}$ can be both positive and negative for packings of circulo-lines. In Sec. III.2, we derive a general expression for the pressure-dependent shear modulus $G$ in terms of derivatives of $U$ and $\phi$ with respect to $\gamma$ at fixed packing fraction and at fixed pressure. Using this expression, we find that the first geometrical family in the $P\rightarrow 0$ limit scales as $G_{f}/G_{0}\sim 1\pm P/P_{0}$ to linear order in $P$, where the sign of $d^{2}\phi/d\gamma^{2}$ determines the sign of the second term in $G_{f}$. In Sec. III.3, we decompose the shear modulus for each low-pressure geometrical family, $G_{f}=G_{a}-G_{na}$, into the affine and non-affine contributions, respectively. We show that the non-affine contribution to the shear modulus of the low-pressure geometrical families is larger for packings of circulo-lines compared to that for disk packings, and that the pressure dependence of $G_{na}$ can cause $G_{f}$ to increase with pressure, which does not occur for packings of spherical particles. In Sec. III.4, we characterize point and jump changes in the contact network during isotropic compression in packings of circulo-lines (interacting via repulsive linear springs) and show that jump changes give rise to discontinuous changes in potential energy, shear stress, and shear modulus, whereas point changes give rise to discontinuous changes only in the shear modulus. We find that point changes are more frequent for packings of circulo-lines (compared to disk packings), but the resulting jumps in the shear modulus are smaller, over the same range of pressure as for disk packings. In Sec. III.5, we discuss the results for the ensemble-averaged shear modulus $\langle G\rangle$ as a function of pressure. We show that, for large aspect ratios ${\cal R}\gtrsim 1.2$, $\langle G\rangle$ scales as a power-law at large pressures with an exponent that is nearly a factor of $2$ larger than that for disk packings. For small aspect ratios, ${\cal R}\lesssim 1.2$, $\langle G\rangle$ versus pressure does not possess a single power-law scaling exponent. We further decompose $\langle G\rangle=\langle G_{f}\rangle+\langle G_{r}\rangle$ into contributions from geometrical families $G_{f}$ and changes in the contact network $G_{r}$. We show that $\langle G_{r}\rangle$ is smaller for packings of circulo-lines compared to that for spherical particles, and thus the increase in the power-law exponent is caused by the fact that the shear modulus can increase in pressure during geometrical families for packings of circulo-lines.

In this article, we studied the structural and mechanical properties of jammed packings of circulo-lines with frictionless, purely repulsive, linear spring interactions. We found several important results for jammed packings of circulo-lines that are different from those for jammed packings of spherical particles. First, we showed that packings of circulo-lines posses geometrical families that can be both concave upward or concave downward in the packing fraction-shear strain ($\phi$-$\gamma$) plane. In contrast, the geometrical families are nearly always concave upward in the $\phi$-$\gamma$ plane, especially at low pressure, for jammed packings of spherical particles. We then derived a stress-dilatancy relation for packings at finite pressure, which allowed us to show that the shear modulus for low-pressure geometrical families obeys $G_{f}=G_{0}+\eta P$ to linear order in pressure, where the sign of $\eta$ is determined by the negative curvature of geometrical families in the $\phi$-$\gamma$ plane. Thus, the shear modulus of low-pressure geometrical families increases with pressure when $d^{2}\phi/d\gamma^{2}<0$ and decreases with pressure when $d^{2}\phi/d\gamma^{2}>0$. The fact that the shear modulus of geometrical families can increase with pressure has a profound effect on the pressure-dependent, ensemble-averaged shear modulus $\langle G(P)\rangle$. In particular, we found that $\langle G(P)\rangle$ for jammed packings of circulo-lines with aspect ratios ${\cal R}\gtrsim 1.2$ displays robust power-law scaling over a wide range of pressure, but the scaling exponent ($\beta\sim 0.8$-$0.9$) is nearly a factor of two larger than that for jammed disk packings ($\beta\sim 0.5$). For smaller aspect ratios, ${\cal R}\lesssim 1.2$, $\langle G(P)\rangle$ does not possess a single power-law scaling exponent over the same range of pressure. To understand the origin of this behavior, we decomposed $\langle G\rangle$ into separate contributions from geometrical families, $\langle G_{f}\rangle$, and from changes in the contact network, $\langle G_{r}\rangle$: $\langle G\rangle/\langle G_{0}\rangle=\langle G_{f}\rangle/\langle G_{0}\rangle+\langle G_{r}\rangle/\langle G_{0}\rangle$, where $\langle G_{0}\rangle$ is the value of $\langle G\rangle$ in the zero-pressure limit. In general, we found that $\langle G_{r}\rangle/\langle G_{0}\rangle$ is larger for disk packings compared to that for packings of circulo-lines, even though the frequency of changes in the contact network is larger for packings of circulo-lines. In contrast, we found that $\langle G_{f}\rangle/\langle G_{0}\rangle$ is much larger for packings of circulo-lines. In fact, $\langle G_{f}\rangle/\langle G_{0}\rangle<0$ for disk packings, whereas it can be positive for packings of circulo-lines in the pressure regime where the power-law scaling exponent is larger than that for disk packings. Thus, the presence of geometrical families with shear moduli that increase with pressure gives rise to important changes in the pressure-dependent mechanical properties for jammed packings of circulo-lines.

These results suggest several promising areas of future research. First, in the present studies, we did not examine in detail how properties of jammed packings of circulo-lines in the ${\cal R}\rightarrow 1$ limit compare to those for jammed disk packings. Two issues arise in the ${\cal R}\rightarrow 1$ limit. 1) Non-circular particles always possess $3$ degrees of freedom per particle, whereas smooth disks possess only $2$ nontrivial degrees of freedom per particle. Thus, in future studies, we will compare the properties of jammed packings of circulo-lines to those for packings of weakly frictional or bumpy particles, which both possess three degrees of freedom per particle [36]. 2) The current force model for circulo-lines, which considers forces between the end and middle sections of pairs of circulo-lines, does not converge to the force model obtained from Eq. 1 in the ${\cal R}\rightarrow 1$ limit for packings with finite pressure. A force law that is a function of the square-root of the area of overlap between pairs of circulo-lines is a more promising model.

In addition, we know that jammed packings of circulo-lines possess concave upward and concave downward geometrical families in the $\phi$-$\gamma$ plane. Do jammed packings of other non-spherical particle shapes also possess concave upward and concave downward geometrical families? Can we find paticle shapes for which jammed packings only possess concave downward geometrical families in the $\phi$-$\gamma$ plane? We have shown in previous studies that the power-law scaling exponent for $\langle G(P)\rangle$ for jammed packings of ellipse-shaped particles is also elevated relative to that for jammed disk packings [25, 26]. Thus, it is likely that jammed packings of ellipse-shaped particles also possess concave downward geometrical families in the $\phi$-$\gamma$ plane. Since the power-law scaling exponent $\beta$ for $\langle G(P)\rangle$ depends on properties of the geometrical families, the frequency of contact network changes, and the size of the discontinuous jumps in $G$ caused by the contact network changes, it seems likely that the power-law scaling exponent $\beta$ will depend sensitively on particle shape. Thus, it will be important to study $\langle G(P)\rangle$ and other mechanical properties for packings of many different particle shapes in both two- and three-dimensions.