Mechanical response of particle packings at jamming onset

Abstract

The jamming onset of particle packings provides an ideal physical state for the comparison of mechanical properties. In this paper, we investigate the effects of size distribution and particle shape on the mechanical response of disordered jammed packings. For packings of bidisperse spheres, both species participate in the force network for the small particle volume fraction X_s ≤ 0.35, boosting the bulk modulus B by nearly 60% relative to the monodisperse case. On the other hand, the high rattler fraction induced by some polydispersity generally reduces B compared to the monodisperse case. Within the family of superellipsoids and spherocylinders studied, packings of nonspherical particles show higher B than that of spheres, with the self-dual ellipsoids achieving the maximum value, where B exhibits robust power-law scaling with pressure. Moreover, although B shows no direct correlation with packing density, it exhibits a robust linear scaling with the coordination number across different particle shapes. This linear relation persists in polydisperse sphere systems when a volume-weighted coordination number is employed to account for the size disparity. These findings advance our understanding of jamming transitions and facilitate rational granular material design.

---

1 Introduction

Athermal particulate materials, such as static packings of granular materials, foams, and emulsion droplets, are ubiquitous in the natural world and industrial applications.^1,2 As one class of nonequilibrium physical systems, granular materials can jam when compressed to sufficiently large packing densities, while the packing structures remain disordered. A series of experimental and numerical simulation results have verified the jammed disordered state with the lowest but robust packing density at jamming onset ϕ_J ∼ 0.64 for spheres in three dimensions (3D).^3–5 At this state, the granular packing develops a mechanically stable structure that can reversibly withstand further external deformation, with nonzero elastic (bulk B and shear G) moduli.⁶

Many previous works have established distinctive structural and mechanical properties of jammed disordered packings. The elastic moduli were found to follow power-law scalings with increasing pressure under isotropic compression, which, in particular, holds for even ordered particulate systems.^7,8 The pressure dependence of B is dominated by the affine response, whereas the pressure dependence of G has significant nonaffine contributions.⁶ It was further shown that, both B and G scale linearly with the mean coordination number z for networks of harmonic central-force springs,⁹ which agrees well with the numerical results observed in sphere packings.⁶ Deviating from monodisperse to polydisperse ones, the sizes of constituent particles can play an influential role in the mechanical properties of the packings.⁸ The vast majority of work so far has concentrated on bidisperse systems by exploring different combinations of control parameters (R,X_s), where R is the small-to-large particle size ratio and X_s is the volume fraction of smaller particles.

Annabattula et al. focused on stress–strain behavior, demonstrating that packing can get “soft” with increasing size disparity and large particle volume fraction.¹⁰ Kumar et al. found that replacing merely 5% volume fraction with smaller particles at fixed size ratio enhances effective bulk modulus by up to 20%.¹¹ Notably, bidisperse systems exhibit a second jamming transition that enriches the jamming phase diagram and introduces unique mechanical properties.¹² It was reported only recently that B initially increases and then decreases with varying X_s, reaching its maximum at X_s ≈ 0.3 for the lowest size ratio considered R = 0.15, which is far from the optimal combination leading to the highest ϕ_J.¹² Compared to bidisperse systems, polydisperse sphere packings exhibit higher shear strength,^13,14 while increasing polydispersity slightly reduces B.^15,16

Particles encountered in nature and engineering practice, however, are normally not spherical but of irregular shapes. Extensive studies have described the dependence of mechanical properties and stress–strain behavior of granular systems on particle shapes. Athanassiadis et al. evidenced that particle shapes can significantly change the effective modulus and the yield stress for a given confinement pressure.¹⁷ Moreover, particular attention has been paid to concave particle packings, which possess both relatively high mechanical stability and low packing densities. Therefore, if carefully designed, concave-shaped particles can serve as promising candidates for light-weight construction materials. Deng et al. systematically investigated the mechanical properties of disordered jammed packings of intersecting spherocylinders.¹⁸ For all the shapes concerned, B is positively correlated with ϕ_J. Interestingly, crystallization can be induced to manipulate and enhance mechanical properties of packings. Compared to disordered jammed packings, triaxial compression tests revealed that these granular crystals are up to 25 times stronger.¹⁹ However, such structures are far beyond the scope of typical granular materials and are not the focus of this study.

Despite studies involving size distribution, shape variations, and coupling effects, inconsistent setup and packing protocols hinder universal characteristics. This underscores the necessity for systematic investigations under an identical packing state, which is the core focus of this work. In this study, we focus on disordered packings at jamming onset and consider both the effects of size distribution and particle shape on B. To ensure consistency and comparability across different systems, all simulations are performed using the linear spring potential. For bidisperse systems with large size disparity, the simultaneous involvement of both large and small particles in the force chains results in an increase in the B by nearly 60%. Compared with monodisperse packings, the high fraction of rattlers induced by particle size polydispersity reduces the system's B.

It is noted that static packings of non-spherical particles at jamming onset typically possess fewer contacts than the isostatic value, i.e., hypostatic.²⁰ How B evolves near the jamming point of non-spherical particle packings, and to what extent hypostaticity affects B, are among the central questions discussed in this work. Various non-spherical particle shapes including superellipsoids, spherocylinders, and their counterparts in two dimensions (2D) are considered. Under fixed confining pressure, packings of superellipsoids and spherocylinders demonstrate that non-spherical shapes substantially increase B relative to spheres, with the self-dual superellipsoid achieving the maximum modulus value. In all cases, B scales linearly with z, and exhibits robust power-law scaling with confining pressure. This linear relation even holds for spherical systems with size distribution when considering a volume-weighted z. Our findings establish z as the determinative parameter governing B, and provide explicit guidance for the design of granular materials.

The remainder of this paper is organized as follows: in Section 2, we describe the particle-size distribution, the superellipsoid model, packing generation protocols, and calculations of the mechanical properties investigated in this work. Then the simulation results are presented and analyzed in Section 3 and concluding remarks are given in Section 4. We also include the details of analytical linear formulation of elastic moduli in Appendix A.

2 System preparation

In this section, we first introduce the polydisperse sphere systems and the superellipsoid model explored in this work, followed by the packing-generation protocols utilized to generate disordered jammed packings. After that, we describe the calculation of mechanical properties.

2.1 Size distributions of spheres

When size distribution is introduced, the packing behavior in polydisperse systems can differ remarkably from that of monodisperse ones. It has been shown that the density increment from 0.64 in polydisperse sphere packings correlates positively with the polydispersity δ,²¹ defined aswhere σ is the particle diameter and 〈·〉 denotes the average over the entire packing system.

We start by considering that the particle diameter follows a log–normal distribution, i.e.,

in which ln σ is distributed according to a Gaussian probability density function (PDF), with the mean 〈σ〉 and the standard deviation . The size polydispersity will vanish at α = 0. δ becomes monotonic with α as . We note that eqn (2) is nearly identical to a normal (Gaussian) distribution for small polydispersities, but has the advantage that it is zero for negative diameters. We also study systems with a uniform size distribution centered around 〈σ〉. We studied polydispersities up to δ = 0.4. For a given δ, particle diameters are sampled uniformly from the interval .

The polydispersity is also sampled from inverse power law (IP, P(σ) ∝ σ⁻³, σ ∈ [σ_min,σ_max]) size distributions. Let k = σ_max/σ_min, then

and we obtain

2.2 The superellipsoid model

Superellipsoid is a general continuous shape family beyond the elementary ones, and is defined aswhere a, b, and c are the lengths of the semi-major axes, and p is the surface shape parameter. a : b : c = α : α^β : 1(a ≥ b ≥ c) gives two independent parameters that control the particle shape: α ≥ 1 represents the aspect ratio and β (0 ≤ β ≤ 1) denotes the oblateness.²² When β = 1/2, the shape corresponds to self-dual ellipsoids. For β = 0 or β = 1, the particle shapes reduce to spheroids characterized by a single aspect ratio w, under which circumstances w = α > 1 for prolate ellipsoids and w = 1/α < 1 for oblate ellipsoids. When a = b = c = 1, the shape degenerates to superballs.²³

We study the mechanical properties of jammed packings of N frictionless superellipsoids that interact via the pairwise, purely repulsive potential²⁴

where ε is the energy scale, η_ij represents the factor that rescales their axes to bring them to exact tangency. Particles i, j overlap when η_ij < 1. The corresponding repulsive force on particle i, arising from particle j, iswhere r⃑_ij is the unit vector pointing from the center of particle j to that of particle i, n̂_ij denotes the unit normal of the tangency point directed toward particle i. The torque is calculated using

[small tau, Greek, vector]_ij = l⃑_ij × f⃑ij, (7)

where l⃑_ij is the moment arm. We set the units of energy and length to be ε and the equivalent volume diameter D_v of constituent particles, respectively. Pressure, shear stress, and elastic moduli are expressed in units of ε/D^d_v, d = 2, 3 is the spatial dimension, while forces are measured in ε/D_v. We note that bidispersity is introduced to suppress ordering for 2D packings, with one-third of particles being 1.4 times larger than the remaining two-thirds,^22,25 and D_v is defined based on the smaller particles. To compare mechanical responses across different shapes, we rescale particle sizes to maintain D_v = 1, ensuring approximately equivalent stiffness.

2.3 Packing-generation protocols

We distinguish between spherical and non-spherical shapes, and use two different protocols to generate their disordered jammed packings.

2.3.1 Discrete element method. In general, the discrete element method (DEM) refers to a molecular dynamics (MD) type approach based on the soft-particle model, typically employing explicit time-integration schemes to solve particle motion; due to its computational simplicity and intuitive implementation, it is widely used for numerical simulations of particulate systems. In DEM, particle trajectories closely resemble reality by following Newton's second law of motion, and the system evolves by (approximately) integrating the equations of motion in time. Moreover, DEM often incorporates damping to dissipate system energy, thereby enabling rapid relaxation toward a stable state.

The present study of spherical particles is limited to frictionless, purely repulsive normal contacts with the linear spring–dashpot model.²⁶ Let m denote the mass of a single particle, then the normal force between contacting particles i and j is given by

f⃑_n = k_nδn̂ − m_effγ_nv⃑_n, (8)

where k_n = 1 is the Hookean spring constant, δ = σ_ij − |r⃑_ij| denotes the particle overlap, m_eff = m_im_j/(m_i + m_j) is the effective mass, and γ_n is the damping coefficient. The unit vector n̂ = r⃑_ij/|r⃑_ij| points from the center of particle j to particle i, and v⃑_n is the relative velocity projected along n̂. Note that we zero out the attractive normal force with limit damping.

LAMMPS is employed in this work to perform DEM simulations using the NPH ensemble.^27,28 An initial dilute configuration, e.g., ϕ = 0.3, of particles with random positions and zero velocities is generated using the random sequential adsorption (RSA) protocol. Thereafter, the packing boundary undergoes isotropic compression to maintain the applied external pressure. We use the velocity-Verlet scheme to integrate particle translational moves, as well as the length of packing boundaries. A constant pressure P_a = 10⁻⁶ (in units of k_n/D, where D is the chosen length scale) is applied together with a pressure damping coefficient P_damp = 2.25τ⁻¹, where defines the characteristic time scale. All particles have material density ρ = 1, the normal damping coefficient is γ_n = 0.5τ⁻¹. The simulation time step is chosen as Δt = 0.02τ.

The onset of jamming in a spherical system with periodic boundary conditions occurs when the number of interparticle contacts N_c first reaches the isostatic value,²⁹N^iso_c = d_fN₀ − d + 1, where d_f = d is the degree of freedom per particle, N₀ is the number of non-rattler particles. Non-rattler particles are those that are considered force-bearing and contribute to the force balance of jammed packings. We define the average number of contacts per particle as z ≡ 2N_c/N₀. Herein, the packing system is deemed force-balanced and we terminate the simulation when the average kinetic energy satisfies U_k/N < 10⁻¹⁶k_nD, and concurrently the average normalized force is below a small threshold, i.e., 〈|f⃑_i|〉/〈|f⃑_ij|〉 < Δ_th, where Δ_th = 10⁻⁴.²⁰ It is noted that the rattler particles are recursively excluded from the above computation.

Notably, the yielding static packings correspond to the particle contact model with

and Θ(·) being the Heaviside step function, so that particles of different sizes exhibit approximately the same effective stiffness in the linear elastic regime of this model. It is also noted that eqn (9) corresponds to the special case of eqn (5) for spherical shapes, differing only by a length-scaling factor. Since stress has units of energy per unit volume, i.e., k_n/D^d−2, the packing system is highly sensitive to volume changes caused by variations in (R,X_s). Therefore, unlike previous studies that fix the size of one component, we set the length unit to 〈σ〉, i.e., D = 〈σ〉, thereby accounting for size distributions.

2.3.2 Energy minimization. For non-spherical particles, appropriate implementation of the contact mechanics is crucial. An athermal quasistatic compression protocol is thus utilized here, and the resulting disordered packings are mechanically stable,^5,30 despite the fact that N_c < N^iso_c. We first initialize an overlap-free, dilute configuration of particles with random positions and orientations, and then successively compress the system using an athermal quasistatic compression protocol with the increment in ϕ_J as Δϕ = 1 × 10⁻⁴.⁶ After each time step, the system is relaxed to the nearest local potential energy minimum using the fast inertial relaxation engine (FIRE) minimization method,³¹ which is extended here to accommodate non-spherical particles.

Specifically, all non-spherical particles (including axisymmetric ones) in three dimensions (3D) are deemed as possessing three rotational degrees of freedom for computational simplicity, i.e., d_f = 6, and the degrees of freedom for particle i are represented as ( for the 2D case, where d_f = 3) for dimensional consistency, where I_i denotes moment of inertia, and α_i, β_i, γ_i are Euler angle components. The energy minimization is terminated when the system has a total net force satisfying , where is the total potential energy. A binary search algorithm is also utilized to push the system to a target pressure P_t = 10⁻⁷ with |P − P_t| < 10⁻⁵P_t. All the results in this work are the averages over 40 independent implementations.

Starting from the packings at jamming onset, we perform minimization of the enthalpy H = U + P₀V, where P₀ is the target pressure, using the FIRE algorithm with a fixed cubic box shape.³² Each of the packings is sampled logarithmically in pressure, spanning from isostatic packings at P₀ = 10⁻⁷ to compressed states with P₀ = 10⁻².

2.4 Structural properties of packings

We calculate the global stress tensor Σ_αβ for each packing sample using the virial expression:where V = L^d is the volume of simulation box, indices α, β label the coordinate axes. Note that eqn (10) reduces to for spherical shapes. We define the shear stress as Σ = −Σ_αβ and the pressure as P = tr(Σ)/d.

For non-spherical particles, we instead employ the Love expression:³³

where p_ij corresponds to the vector from the center of particle i to the contact point C_ij with particle j. We also note that the Love expression reduces to the virial expression for spherical particles.

The fabric tensor is a key descriptor of the orientational distribution of interparticle contacts, and is defined as³⁴

where V_p is the volume of particle p and n̂^k is the unit normal from the center of p to its k-th contact point. Note that can be viewed as the volume-weighted average of z.³⁴ We term tr(F) as structural isotropy F_v.

To quantify structural anisotropy, we also use the fabric deviator:³⁵

where λ₁ ≥ λ₂ ≥ λ₃ are the eigenvalues of F. A larger λ_d indicates a more anisotropic contact network.

2.5 Calculation of mechanical properties

The linear elastic response of the isotropic systems can be understood by characterizing their ability to resist global affine deformation. Specifically, B is for volume-changing bulk deformation ε_B, where ε represents the strain tensor. An important feature of the mechanical response of jammed particle packings is the dependence of elastic modulus on P under isotropic compression. The elastic modulus can be decomposed into the affine and nonaffine contributions, i.e., B = B_a − B_na. The affine contribution considers the linear response of the jammed packing to an ideal simple shear deformation without relaxation, whereas the nonaffine contribution includes particle motion from minimization of the total potential energy. Instead of explicitly applying a deformation field, we calculate the elastic moduli by imposing virtual deformations and calculating the derivatives of energy, of which the details are provided in Appendix A.

For non-spherical particles, however, we numerically compute elastic moduli via infinitesimal deformations on packings instead due to complications of particle interactions. To be specific, volumetric strain δL/L = 10⁻⁶ is applied after obtaining the disordered jammed packings, and we clarify that the contact networks in packings do not change during the deformation. Post-deformation energy minimization captures non-affine displacements, yielding

3 Results

We present our results and discussion in three subsections. In Section 3.1, we investigate disordered jammed packings of spheres, including bidisperse and three representative polydisperse cases, and systematically analyze elastic moduli near point J. In Sections 3.2 and 3.3, we turn to monodisperse disordered jammed packings of two classes of convex particles: superellipsoids (including their 2D degenerations) and 2D spherocylinders, thereby covering representative nonspherical shapes via parameter tuning. We also show the pressure-dependent scaling of B.

3.1 Polydisperse packings of spheres

As the simplest polydisperse case, bidisperse sphere systems have received extensive attention and study. We consider a system size of N = 6000 particles, and perform DEM simulation under fixed pressure P = 10⁻⁶. The number of large particles decreases rapidly as X_s approaches unity. Hence, for computational efficiency, the discussion is limited to R ≥ 0.1, and results are averaged over 10 independent implementations.

In Fig. 1, we show the global bulk modulus B of bidisperse spheres analytically as the function of X_s, which differs from local moduli discussed in ref. 12. For small R, as X_s increases from zero, B initially decreases, and then increases to a maximum at X_s ≈ 0.4. Using B ≈ 0.26 for monodisperse cases as a reference, the regime 0 < X_s < 0.35 (region I) for bidisperse spheres corresponds to a regime of reduced B, while B rises for X_s ≥ 0.35 (region II). By contrast, for a larger ratio R = 0.85, B remains unchanged as X_s varies.

Fig. 1 Bulk modulus B of bidisperse sphere packings versus the volume fraction X_s for different size ratios R.

Furthermore, we compute the fraction of rattlers N_r/N with N_r = N − N₀, as shown in Fig. 2(a). For small R, nearly 99% of constituent particles become rattlers within region I. In region II, N_r/N drops significantly below 10%, and approaches around 1.5% as X_s increases, with the transition point at approximately X_s ≈ 0.3. Fig. 2(b) and (c) show the disordered jammed packing structures at R = 0.1 in region I and region II, respectively. Clearly, only large particles form the force chain in region I, while both large and small particles participate in the force chain, yielding an increase in B. As R decreases, region II exhibits a lower N_r/N, which explains its higher B. Additionally, for small values of R and X_s, N_r/N undergoes a transition when the force network shifts from being dominated by large–large contacts to one in which both large and small particles participate, resembling the behavior observed for B.

Fig. 2 (a) The fraction of rattlers N_r/N in disordered jammed packings of bidisperse spheres versus the volume fraction X_s at different size ratios R. The dashed line corresponds to the case for monodisperse spheres. Disordered jammed packing structures of bidisperse spheres with (R,X_s) = (0.1,0.1) and (0.1,0.4) are shown in (b) and (c), respectively. Each particle is colored according to its local contact number z_i (see the color bar), while rattlers are marked in white.

In Fig. 3, we show the structural anisotropy λ_d of different bidisperse systems. Regions I and II are again distinguished in the curves: in region I, λ_d decreases significantly, while in region II it remains nearly constant. We also examine the variation of B with F_v in region II, as shown in Fig. 4. Overall, the results indicate a clear linear correlation between B and F_v across different R. As for large R = 0.85, both B and F_v remain nearly the same. We note that there is a consistent downturn trend in B for large F_v, with the reduced values corresponding to X_s = 0.4, i.e., the approximate turning point of region II. Owing to finite size effects, packings at X_s = 0.4 are more prone to structural heterogeneity arising from the uneven distribution of larger particles, thus resulting in a lower B.

Fig. 3 Structural anisotropy λ_d for disordered jammed packings of bidisperse spheres versus volume fraction X_s at different diameter ratios R.

Fig. 4 Bulk modulus B in disordered jammed packings of bidisperse spheres versus structural isotropy F_v at different diameter ratios R.

Fig. 5 shows the relationship between B and polydispersity δ for three distributions investigated in this work. In 3D, B decreases monotonically with increasing δ, while the trend is reversed in 2D. Among the three distributions considered, the log–normal distribution exhibits the lowest minimum-to-maximum size ratio, which in turn results in the lowest B in 3D. Note that δ = 0 corresponds to monodisperse packings, which are highly ordered in 2D. The introduction of size polydispersity thus takes on the role of suppressing crystallization. Owing to the absence of a sudden increase in the N_r/N ratio as observed in 3D systems, B in 2D packings increases in a more gradual manner. Further, for 2D polydisperse systems, we exclude packings with high orientational order ψ₆ > 0.1 with ψ₆ being defined as^36,37

where N_j is the number of Voronoi neighbors of particle j, and θ_jk is the angle between an arbitrary fixed axis and the vector connecting centers of particles j and k. As shown in the inset of Fig. 5(a), B in 2D polydisperse systems nearly collapses as a linear function of F_v regardless of the distribution types, similar to the cases of bidisperse sphere systems.

Fig. 5 Bulk modulus B for disordered jammed packings in (a) two dimensions and (b) three dimensions versus polydispersity δ for the three distribution types (uniform, log–normal, inverse power). Inset shows the relationship between B and structural isotropy F_v in two dimensions.

3.2 Disordered packings of superellipsoids

By tuning three shape parameters (p,α,β) in the superellipsoid model (eqn (4)), we investigate over 200 distinct shapes. Fig. 6 displays disordered jammed packings of nine representative superellipsoids. Fig. 7 and 8 present the relationships between control parameters and both ϕ_J and B for uniaxial ellipsoids, superballs, and their 2D counterparts. We limit our results to p ≤ 2.0 due to computational efficiency. Consistent shape-dependent trends for ϕ_J and B are observed across different spatial dimensions. As particle shapes deviate from spherical ones, B initially increases and then decreases, with peak positions differing from those of ϕ_J. Prolate ellipsoids and ellipses achieve maximal B at w ≈ 1.8, while oblate ones peak at w ≈ 0.5. B for superballs and superdisks reach maxima at p ≈ 0.7 and 1.6, whereas ϕ_J monotonically increases with |p − 1| when p ≤ 2.0. Notably, oblate ellipsoids show superior B with respect to prolate ones despite comparable ϕ_J, exhibiting 130% higher B than monodisperse spheres.

Fig. 6 Disordered jammed packings of representative superellipsoids, with control parameters (α,β,p) = (a) oblate ellipsoid (3,1,1), (b) prolate ellipsoid (3,0,1), (c) self-dual ellipsoid (1.7,0.5,1), (d) generic ellipsoid (2,0.25,1), (e) generic ellipsoid (2,0.75,1), (f) superball (1,1,2), (g)–(i) three superellipsoids: (2,1,0.75), (2,0.5,1.1), and (1.5,0.5,1.5).

Fig. 7 Disordered jammed packings of bidisperse ellipses and monodisperse ellipsoids. (a) Packing density ϕ_J and (b) bulk modulus B versus aspect ratio w.

Fig. 8 Disordered jammed packings of bidisperse superdisks and monodisperse superballs. (a) Packing density ϕ_J and (b) bulk modulus B versus surface shape parameter p.

Fig. 9(a) presents the variation of B as a function of control parameters for disordered jammed packings of ellipsoids. In each case, B increases from the value for monodisperse spheres, attains its maximum ∼0.72 for the self-dual ellipsoid (α ≈ 2.25), and then decreases for larger α. By contrast, ϕ_J reaches its peak (ϕ_J ≈ 0.74) at α ≈ 1.5 and subsequently declines, indicating that ϕ_J does not necessarily contribute to B. To decouple the influence of ϕ_J on B, we define the dimensionless modulus B̃ = B/ϕ_J, which measures the stiffness per unit packing density. As shown in Fig. 9(b), B̃ rises monotonically with α until it plateaus near unity as α increases. Notably, the results also confirm the exceptional mechanical efficiency of self-dual ellipsoids at α ≈ 2.25.

Fig. 9 (a) Bulk modulus B and (b) dimensionless modulus B̃ = B/ϕ_J for disordered jammed packings of ellipsoids versus aspect ratio α with varying oblateness β.

Interparticle contact forces are the microscopic origin of macroscopic stress transmission. Fig. 10(a) plots the PDF of the normalized normal force f_n/〈f_n〉 for the self-dual ellipsoids, with the majority of forces distributed between 0.5 〈f_n〉 and 2.0 〈f_n〉, peaking near the mean. Above-average forces decay exponentially, while sub-mean forces follow a power law:

Narrowing PDFs and increasing weak-force exponents at α ≈ 2.25 correlate with force homogenization and maximal B. Insets show right-shifted PDF peaks indicating stress concentration. Fig. 10(b)–(e) display packing structures with rising z versus α, reflecting z–B correlation.

Fig. 10 (a) Probability density function (PDF) of normalized contact force f_n/〈f_n〉 for self-dual ellipsoids at different aspect ratio α. (b)–(e) Disordered jammed packing structures for α = 1.1, 1.5, 2.25, and 3.0, respectively, colored by local coordination number z_i (see color bar). Rattlers (if present) are marked in white.

By tuning the shape parameter p, we obtain more 2D superellipses and 3D superellipsoids. Compared to p = 1 cases, these particle shapes yield higher ϕ_J, partially attributed to increased cubic ordering in their disordered jammed configurations.²⁰z of these hypostatic packings span the range from 2d to 2d_f. Fig. 11 demonstrates the relationship between B and z across all disordered jammed packings. Both 2D and 3D systems show monotonic increases in B with z, following approximately linear trends. Notably, modifying shape parameters from self-dual ellipsoids to cuboid-like superellipsoids (p = 1.5 and 2) reduces both z and corresponding B, confirming self-dual ellipsoids as the optimal shape for maximal B within the superellipsoid family. Similarly, ellipses achieve peak B among the 2D superellipses considered. Defining reduced bulk modulus B/B₀ normalized by spherical/disk systems, we observe the proportionality coefficient in the linear relationship approximates the inverse of spatial dimensionality d, i.e.,

where B/B₀ > 1 universally. Self-dual ellipsoids achieve maximal B among the superellipsoids, paralleling ellipses in 2D.

Fig. 11 Reduced bulk modulus B/B₀ for disordered jammed packings of 2D superellipses and 3D superellipsoids versus coordination number z.

We further examine the power-law scaling of elastic moduli with pressure P for superellipsoids, and generate 50 compressed configurations per initial packing. Fig. 12 shows representative pressure-dependent elastic moduli for superellipsoid packings. Within 10⁻⁷ < P < 10⁻⁴, B exhibits a plateau regime followed by progressive stiffening. Defining the plateau modulus B_c, the inset demonstrates consistent power-law scaling B–B_c ∼ P^β, and we further verify that β ≈ 1.0 is universal across different shapes among the family of superellipsoids.

Fig. 12 Pressure P dependence of bulk modulus B for disordered jammed packings of self-dual ellipsoids. Inset shows the relationship between B–B_c and P with B_c being the plateau value of B.

3.3 Disordered packings of spherocylinders

Despite similar elongation effects, spherocylinders differ fundamentally from ellipsoids. Following the discussion of ref. 38, contacts between different spherocylinders can be classified into three categories: end–end, end–middle, and middle–middle. Linear spring potential as shown in eqn (5) is adopted for the three types of contacts, and it is noted that middle–middle contacts are modeled as combinations of two end–middle contacts to ensure continuity in both U and f. Also, middle–middle contacts contribute two contacts to the computation of z.

Using FIRE minimization, we generate jammed packings of 2D bidisperse spherocylinders (diameter σ) under P = 10⁻⁷. Fig. 13(a) displays B versus aspect ratio w for 2D spherocylinder packings. B exhibits a non-monotonic trend with w, peaking at w = 1.4, while ϕ_J reaches its maximum near w ≈ 2. Notably, spherocylinders demonstrate lower B values compared to ellipses. The inset reveals the relationship between reduced bulk modulus B/B₀ and z: linear scaling (as described by eqn (17)) holds for small z (w ≤ 1.4), transitioning to monotonic decline with increasing z at w > 1.4. This dichotomy stems from contact mechanics: middle–middle contacts are virtually absent below w ≤ 1.4, making z equivalent to average contact number. Beyond this threshold, middle–middle contacts progressively dominate. Essentially, near-spherical spherocylinders (w ≤ 1.4) exhibit superellipsoid-like contact patterns governed by eqn (17), whereas elongated spherocylinders (w > 1.4) develop static determinacy through end–end contacts counted as dual constraints.³⁹Fig. 13(b) and (c) contrast local structures at w = 1.2 (sphere-like contact network) and w = 2.0 (end-contact dominated configuration), highlighting fundamental differences in force transmission mechanisms.

Fig. 13 (a) Bulk modulus B for disordered jammed packings of 2D spherocylinders versus aspect ratio w. Inset shows the relationship between reduced bulk modulus B/B₀ and coordination number z. Local views of disordered jammed packings for w = 1.2 and w = 2.0 are shown in (b) and (c), respectively.

4 Conclusions

We systematically investigate disordered jammed packings and their mechanical response for monodisperse, bidisperse, and polydisperse sphere systems. In bidisperse sphere packings, we introduce size variation and identify two distinct jamming modes: one dominated by contacts among large particles, and the other involving both large and small particles in the force chain. As the small-particle volume fraction X_s and the diameter ratio R vary, the elastic moduli display a complex, non-monotonic dependence, reflecting these two jamming scenarios. When X_s < 0.35, small spheres fill the voids among large spheres, yielding a high rattler fraction and a softened bulk modulus B. In contrast, at higher X_s, both species contribute to the load-bearing network, and B increases markedly. This finding has implications for understanding the mechanical behavior in bidisperse and more general polydisperse granular systems. At the same time, size disparity in the bidisperse system does not increase structural anisotropy but rather enhances structural isotropy F_v, and B is found to vary linearly with F_v, underscoring the critical role of coordination-like parameters in controlling the modulus. By tuning (X_s,R), one can effectively adjust B, achieving nearly a 60% increase relative to the monodisperse reference. For polydisperse spheres, we consider uniform, inverse power, and log–normal size distributions. In three dimensions, the global bulk modulus B decreases monotonically as the polydispersity δ increases; in two dimensions, the trend is reversed and B increases with δ. This contrast highlights the important influence of the rattler fraction and suppression of crystallization under size polydispersity. The linearity between B and F_v is also observed for 2D polydisperse systems with the prerequisite of structural disorder. These results provide strong support for understanding how microscopic contact mechanisms in different packing structures affect the macroscopic elastic response.

We further systematically examine the effect of non-spherical particle shapes on elastic moduli in disordered jammed packings. We consider representative families including 3D superellipsoids, 2D superellipses, and 2D spherocylinders. Numerical results indicate that as particle shape gradually deviates from perfect sphericity, the bulk modulus B first increases and then decreases, with the peak of B occurring at a shape parameter different from the peak of the jamming packing fraction ϕ_J, showing that B and ϕ_J are not in one-to-one correspondence. Using the monodisperse sphere (or bidisperse disk in 2D) modulus B₀ as a baseline, we find B/B₀ > 1 throughout the examined shape range, indicating superior compressive stiffness of non-spherical particles compared to spheres. Among 3D superellipsoids, the self-dual shape achieves the highest B, with a peak increase of nearly 130% at the optimal shape parameter α = 2.25, while in 2D superellipses the optimal shape is an ellipse. Analysis of force-chain networks reveals that the increase in B is associated with both enhanced stress concentration and a rightward shift of the peak in the stress distribution.

We further uncover that in near-spherical disordered jammed packings, the ratio B/B₀ depends approximately linearly on the average coordination number z in a manner independent of particle shape, with the slope equal to the reciprocal of the spatial dimension. This conclusion corroborates the optimality of the self-dual superellipsoid and aligns with the linear B–F_v relation found in bidisperse sphere systems, jointly demonstrating that z is the decisive parameter for B. Upon further isotropic compression of the obtained superellipsoid jammed structures, we observe that B exhibits robust power-law scaling with P. In the high-pressure regime, the exponent for B approaches unity, independent of particle shape and spatial dimension. For two-dimensional spherocylinders near the spherical limit, similar behavior to superellipses is observed, such as the linear B/B₀–z relation, but due to differences in contact geometry, full unification with superellipsoid results is not achieved and warrants further investigation.

This work systematically elucidates the mechanical characteristics of disordered jammed packings, revealing how size distribution, shape parameters, z, and contact force distributions jointly influence the elastic moduli. These findings provide a solid foundation for understanding microscale contact mechanisms and the resulting macroscopic mechanical behavior in non-spherical particle systems can help guide the design of granular materials via optimizing z. It will be interesting to further discuss the effect of different scaling exponents of the interaction, e.g., the hertz-like model.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article are available at Science Data Bank at https://doi.org/10.57760/sciencedb.28480.

Appendix A: linear response

In this Appendix, we calculate elastic modulus for jammed disordered packings of spherical particles. Formally, one can construct the dN × dN dynamical matrix H_iα,jβ (the indices i and j referring to the particles and α,β = x,y representing the Cartesian coordinates) by taking the second derivative of the total potential energy: H_iα,jβ = ∂²U/∂r_iα∂r_jβ, which can be further rewritten as⁴⁰andwhere δ_αβ is the Kronecker delta.

Following the discussions by Maloney and Lemaitre,⁴¹ the affine modulus M_a is formulated as the second derivative of U with respect to the homogeneous affine strain, i.e.,

In particular, for volumetric deformations B_a requires the following correction:where P = −∂U/∂V.

In disordered jammed systems, an affine deformation typically induces force imbalances, which are usually accompanied by a secondary, non-affine response. Using the dynamical matrix, the additional non-affine displacement Δ[small xi, Greek, vector]_M can be obtained by solving the linear response equation:

where can be interpreted as the external forces resulting from an elementary affine deformation ε_M. Accordingly, we defineDue to the presence of non-affine relaxation, the elastic modulus of a disordered jammed packing is always smaller than the purely affine term. As deformation accumulates, the non-affine response becomes more pronounced, the system approaches instability, and the modulus decreases accordingly. By summing affine and nonaffine contributions, B can be written asNote that this method does not explicitly impose a physical strain field on the system; instead, we apply a virtual affine deformation and compute energy derivatives to determine the elastic moduli.

Acknowledgements

We thank Weiwei Jin and Shiyun Zhang for many insightful discussions and useful advice on the manuscript. We also thank Prof. Alessio Zaccone for helpful information. This work was supported by the National Natural Science Foundation of China (grant no. 11972047) and the High-Performance Computing Platform of Peking University.

Notes and references

1. M. van Hecke, J. Phys.: Condens. Matter, 2009, 22, 033101.
2. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347–369.
3. G. D. Scott and D. M. Kilgour, J. Phys. D: Appl. Phys., 1969, 2, 863.
4. S. Torquato, T. M. Truskett and P. G. Debenedetti, Phys. Rev. Lett., 2000, 84, 2064.
5. C. S. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2002, 88, 075507.
6. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
7. C. P. Goodrich, A. J. Liu and S. R. Nagel, Nat. Phys., 2014, 10, 578–581.
8. H. Mizuno, K. Saitoh and L. E. Silbert, Phys. Rev. Mater., 2020, 4, 115602.
9. A. Zaccone and E. Scossa-Romano, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 184205.
10. R. K. Annabattula, Y. Gan and M. Kamlah, Fusion Eng. Des., 2012, 87, 853–858.
11. N. Kumar, V. Magnanimo, M. Ramaioli and S. Luding, Powder Technol., 2016, 293, 94–112.
12. J. C. Petit, N. Kumar, S. Luding and M. Sperl, Phys. Rev. E, 2022, 106, 054903.
13. M. Wackenhut, S. McNamara and H. Herrmann, Eur. Phys. J. E: Soft Matter Biol. Phys., 2005, 17, 237–246.
14. C. Voivret, F. Radjai, J.-Y. Delenne and M. S. El Youssoufi, Phys. Rev. Lett., 2009, 102, 178001.
15. J. Wiacek and M. Molenda, Granular Matter, 2018, 20, 17.
16. V. Magnanimo, L. La Ragione, J. T. Jenkins, P. Wang and H. A. Makse, Europhys. Lett., 2008, 81, 34006.
17. A. G. Athanassiadis, M. Z. Miskin, P. Kaplan, N. Rodenberg, S. H. Lee, J. Merritt, E. Brown, J. Amend, H. Lipson and H. M. Jaeger, Soft Matter, 2014, 10, 48–59.
18. W. Deng, L. Liu, Y. Yuan and S. Li, Powder Technol., 2021, 383, 443–453.
19. A. N. Karuriya and F. Barthelat, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2215508120.
20. Y. Yuan, K. VanderWerf, M. D. Shattuck and C. S. O'Hern, Soft Matter, 2019, 15, 9751–9761.
21. R. S. Farr and R. D. Groot, J. Chem. Phys., 2009, 131, 244104.
22. A. Donev, R. Connelly, F. H. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 051304.
23. Y. Jiao, F. H. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041304.
24. Y. Yuan, Y. Jiao, Y. Wang and S. Li, Phys. Rev. Res., 2021, 3, 033084.
25. Z. Huang, W. Deng, Y. Yuan, L. Liu, Y. Wang and S. Li, Powder Technol., 2022, 396, 565–577.
26. L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, D. Levine and S. J. Plimpton, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 051302.
27. A. P. Santos, D. S. Bolintineanu, G. S. Grest, J. B. Lechman, S. J. Plimpton, I. Srivastava and L. E. Silbert, Phys. Rev. E, 2020, 102, 032903.
28. J. M. Monti, J. T. Clemmer, I. Srivastava, L. E. Silbert, G. S. Grest and J. B. Lechman, Phys. Rev. E, 2022, 106, 034901.
29. A. V. Tkachenko and T. A. Witten, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 687.
30. C. F. Schreck, M. Mailman, B. Chakraborty and C. S. O'Hern, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 061305.
31. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
32. S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. Van Hecke, Phys. Rev. Lett., 2012, 109, 095703.
33. C. F. Schreck, N. Xu and C. S. O'Hern, Soft Matter, 2010, 6, 2960–2969.
34. N. Kumar, S. Luding and V. Magnanimo, Acta Mech., 2014, 225, 2319–2343.
35. D. Barreto, C. O'Sullivan and L. Zdravkovic, AIP Conf. Proc., 2009, 1145, 181–184.
36. D. R. Nelson and B. Halperin, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 19, 2457.
37. Y. Yuan, W. Jin, L. Liu and S. Li, Phys. A, 2016, 461, 747–755.
38. J. Zhang, K. VanderWerf, C. Li, S. Zhang, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2021, 104, 014901.
39. K. VanderWerf, W. Jin, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2018, 97, 012909.
40. A. Tanguy, J. Wittmer, F. Leonforte and J.-L. Barrat, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 174205.
41. C. E. Maloney and A. Lemaitre, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 016118.

---