Crystallization in load-controlled shearing flows of monosized spheres

Abstract

Identical, inelastic spheres crystallize when sheared between two parallel, bumpy planes under a constant load larger than a minimum value. We investigate the effect of the inter-particle friction coefficient of the sheared particles on the flow dynamics and the crystallization process with discrete element simulations. If the imposed load is about the minimum value to observe crystallization in frictionless spheres, adding small friction to the granular assembly results in a shear band adjacent to one of the planes and one crystallized region, where a plug flow is observed. The ordered particles are arranged in both face-centered cubic and hexagonal-closed packed phases. The particles in the shear band are in between the crystalline state and the fluid state, but the latter is never reached, which results in a large shear resistance. As the particle friction increases, the shear band disappears, and the ordering in the core region is destroyed. A significant portion of the particles are in a fluid state with a zero shear rate, leading to a substantial and unexpected reduction in the shear resistance with respect to the frictionless case. If the imposed load is increased well above the minimum from the onset of crystallization, we observe the formation of one shear band in the core, where the particles are again between the crystalline state and the fluid state, surrounded by two crystallized regions near the boundaries, in which most of the particles are in the face-centered cubic phase and translate as a rigid body with the boundaries themselves. In this case, the macroscopic shear resistance is independent of the particle friction.

---

1. Introduction

Understanding granular flows is crucial for a large number of industrial and geophysical applications. Granular materials exhibit complex behaviour,¹ and they even undergo a phase transition when subjected to external energy, with their structure changing from a disordered state to an ordered state.² Many studies have investigated the effect of the driving force on the crystallization properties of granular matter, like vibration, shearing or their combination. However, their results are controversial: in experiments on continuous shear, the particles were observed to always exhibit crystallization,^3,4 whereas there was no detectable ordering when vibration was applied.^5,6 Experiments on the combined effect of vibration and shearing, on the other hand, have shown opposite results, with vibration inducing and shearing destroying crystallization.⁷ Ciamarra et al.⁸ have observed that both vibration and shear prevent the crystallization of the system provided they are strong enough. Other works on vibration have emphasized the importance of frequency and amplitudes in the crystallization process and the associated structures.^9–16

The presence of crystalline structures can significantly affect the rheology,^2,17 but existing continuum approaches for dense particle flows do not take into account topological ordering.^18,19 For this reason, investigating the ordering at the particle level is crucial. Panaitescu et al.²⁰ studied the nucleation process of shear-induced crystallization in granular sphere packings and revealed the shape and size of the microcrystals. The size and structure of the nuclei were studied in detail by laser tomography experiments.²¹ Duff and Lacks proposed a mechanism based on a fold catastrophe of the free-energy landscape for shear-induced crystallization of jammed states.²² Focusing on the influence of strain and shear rate on the crystallization, Mokshin and Barrat²³ reported two opposite effects of shearing: aiding the formation of small crystallites and suppressing the creation of large clusters. Mesoscopic evolution and properties of the crystal structures in dense granular flow under continuous shear have been explored through experiments and simulations by Bai et al.¹⁸

The geometry of the boundaries significantly affects the crystallization process. Silbert et al.²⁴ systematically investigated the effect of the boundary roughness on the order–disorder transition in inclined granular flows by discrete element simulations with periodic boundary conditions. The ordered states showed layers of particles, arranged in a hexagonal close-packing, parallel to the flow direction and sliding over each other. This results in a lower shear resistance compared to the disordered states. Kumaran and Maheswari²⁵ and Kumaran and Bharathiyar²⁶ confirmed the observation and found that there is a minimum base roughness above which the order–disorder transition takes place. The characteristics of an inclined flow over a base made bumpy with a layer of glued spheres arranged in a triangular lattice, which shows intermittency between the ordered and disordered states, were studied by Yang et al.²⁷

Discrete element (DEM) numerical simulations of pressure-controlled flows of identical, frictionless spheres, sheared between bumpy planes made of regularly arranged frictionless particles, in the absence of gravity,²⁸ have shown that there exists a minimum load above which the spheres crystallize. Here, we want to extend this work and investigate the role played by the inter-particle friction coefficient and the applied load on the crystallization process and the resulting micro-structures.

The paper is organized as follows: in Section 2 we describe the discrete numerical simulations, the associated parameters and the quantities that we have employed to characterize the micro-structures of the flow. In Section 3, we show and comment on the results of the simulations. Final remarks and plans for future works are provided in Section 4.

2. Simulation method

We have performed DEM simulations of granular shear flows with the molecular dynamics platform LAMMPS.† ²⁹ In the simulations, N identical soft spheres of diameter d and mass density ρ_p are initially placed at random in a rectangular box of length L_x = 20d, width L_y ≈ 10d, and height L_z ≈ 20d and then sheared between two parallel planes moving relative to each other with a constant velocity 2V, in the absence of gravity. The x, y and z directions are the flow, gradient and vorticity directions, respectively. The planes are constructed by glueing one layer of particles having the same diameter as the flowing particles in a hexagonal closed-packed arrangement. Each plane consists of N_w spheres. Both planes experience a constant compressing load F_z along the z-axis, and can freely move vertically as rigid bodies. We apply periodic boundary conditions along the x and y directions. A snapshot of the shear cell with the associated frame of reference is shown in Fig. 1.

Fig. 1 A snapshot of the DEM simulation: the flowing particles (coloured in yellow) are sheared by moving the two planes in the opposite direction along the x-axis at a constant relative velocity of 2V.

We employ the Hookean contact force model. In this model, the two contacting particles i and j (at positions [r_i, r_j], with diameters [d_i, d_j], masses [m_i, m_j], translational velocities [v_i, v_j], and angular velocities [ω_i, ω_j]) experience a relative normal compression with overlap δ = d_i + d_j − r_ij, where r_ij = r_i − r_j and r_ij = |r_ij|. The force on particle i due to its interaction with particle j, F_ij, is the sum of normal and tangential contributions: F_ij = Fⁿ_ij + F^t_ij, given as:

Fⁿ_ij = K_nδn_ij − γ_nv_n, (1)

F^t_ij = −K_tΔs_t − γ_tv_t, (2)

where n_ij = r_ij/r_ij; v_n and v_t are the normal and the tangential relative velocities between the particles i and j, respectively; K_n and K_t are the elastic constants, whereas γ_n and γ_t are the visco-elastic damping constants, respectively. The values of γ_n, γ_t and K_t are determined as follows:andwhere e_n is the normal coefficient of restitution, i.e., the negative of the ratio of the normal relative velocities after and before the collision; e_t is the tangential coefficient of restitution, i.e., the negative of the ratio of the tangential relative velocities after and before the collision; and m_ij = m_im_j/(m_i + m_j) is the effective mass.

The collision time can be found as

The quantity Δs_t in eqn (2) denotes the elastic tangential displacement accumulated during the entire duration of the contact between spheres. The magnitude of tangential displacement Δs_t is truncated in order to satisfy a local Coulomb yield criterion: |F^t_ij| < μ_ij|Fⁿ_ij|, where , and μ_i and μ_j are the friction coefficient of particles i and j, respectively.

The total force and torque on particle i are obtained as

where the summation is carried out over all j particles in contact with particle i. The translational and rotational motion of particle i are then determined by Newton's second law of motion:where I_i is the moment of inertia of the particle i. The position and velocity of the particle i are updated by a numerical integration at every time step. The contact laws for interactions between a flowing particle and a boundary particle are the same as those between two flowing particles. Note that since the boundary particles are subjected to an external force in addition to the contact forces, we use the LAMMPS “fix aveforce” command to determine their motion. The total force on the plane is measured by summing forces over the boundary particles and the imposed load F_z:which then is averaged over the N_w boundary particles so that each of them experiences the same force (F_w/N_w), resulting in a rigid body motion of the plane over time.

All the quantities are made dimensionless using d, ρ_p and V, so that distances, times, velocities, forces, elastic and viscoelastic constants are measured in units of d, d/V, V, ρ_pd²V², ρ_pdV², and ρ_pd²V, respectively. All the particles have the same mass, m = π/6. The boundary particles are frictionless, i.e., their friction coefficient is μ_w = 0. We have employed N = 3150 flowing particles and N_w = 240 boundary particles per plane, with e_n = e_t = 0.4 and K_n = 2 × 10⁵. The integration time step δt for the numerical integration is fixed and equal to t_c/50, that is δt ∼ 7.5 × 10⁻⁵. The saving time step of the simulations is set to 100 integration time steps, so that measurements are recorded at time intervals of 0.0075.

We have chosen seven values of the friction coefficient of the moving particles, μ_p, in the range between 0 and 0.5 for two different applied loads, F_z = 2.4 and F_z = 24 000. The lower value of the load is about the minimum value required to observe crystallization in imposed-pressure shearing flows of frictionless spheres when e_n = 0.4.²⁸ The simulations require running about 10⁹ time-steps to reach a steady state for F_z = 2.4 and about half of them for F_z = 24 000.

We have measured the temporal evolution of the average solid volume fraction, [small nu, Greek, macron], defined as the fraction of the total volume occupied by the flowing particles, i.e. [small nu, Greek, macron] = Nπ/(6V_t), where π/6 is the volume of a flowing sphere and V_t = L_xL_yH is the volume comprised between the two moving planes. Here, H is the flow height, that is the gap between the edges of the top and bottom boundary particles. We have also measured the temporal evolution of the velocity of the center of mass of the moving particles in the x-direction, , and the average kinetic energy per particle, , where v_i,x is the x-component of the velocity, v_i, of particle i and ‖‖ denotes the Euclidean norm of a vector. We have employed these quantities to assess whether the simulation has reached a steady state (average kinetic energy per particle independent of time), whether the steady state is symmetric (the velocity of the center of mass of the particles is zero) and to estimate whether crystallization has likely occurred (average solid volume fraction larger than the random close packing, 0.64).

Once the steady state is attained, we have measured the local profiles along the z-direction of the solid volume fraction, ν, the x-component of the flow velocity, v_x, and the granular temperature, T, i.e., one third of the mean square of the particle velocity fluctuations around the mean. To do that, we have used the post-processing tool “fstatistics” provided with the open-source DEM code MercuryDPM‡.^30,31 For all the simulations, we have divided the domain comprised between the two planes in 24 equal-sized bins parallel to the x–y plane, and we have then coarse-grained within each bin to obtain local profiles of the continuum variables. Finally we have performed time averaging, over at least 10⁶ saving time steps.

To quantify the crystalline structures in the flow, we could employ bond-orientational order parameters, which were introduced by associating with every bond joining a particle and its neighbours a set of spherical harmonics:³²

where Y^m_l is spherical harmonics and θ_ij and ϕ_ij denote the polar and azimuthal angles which define the orientation of the bond vector pointing from particle i to its neighbour particle j; and N_b(i) contains the set of nearest neighbour indices for particle i. However, the order parameter q_lm is not rotationally invariant. Eslami et al.³³ have proposed the following expression for a rotationally invariant and local order parameter:where is the complex conjugate of and is defined as follows:

Then, the local order parameter can be averaged over the first coordination shell neighbours of particle i³³ to permit a more accurate determination of the local crystalline structure:

The time-average of the local order parameter in the steady state, indicated as , permits to further narrow its distribution. Eslami et al.³³ have shown that the combination of and is sufficient to discriminate between a fluid (disordered) and a solid (crystalline) phase. It can also distinguish within the crystalline phase between body-centered cubic (bcc), face-centered cubic (fcc) and hexagonal close-packed (hcp) crystal structures. For calculating the local order parameters, we choose a cutoff distance of 1.4d to determine N_b(i) ³³ and exclude the boundary particles from the neighbour sets.

Finally, we have measured the macroscopic friction coefficient μ in the steady state as the time-averaged value of the ratio of the tangential force to the normal force exerted on the moving planes. This quantity is representative of the shear resistance of the granular flow and has significant practical implications.

3. Results

We first show and comment on how the inter-particle friction coefficient affects the crystallization process and the shear resistance under continuous shearing and a constant normal load F_z = 2.4. Then, we report how these results change upon increasing the load by four orders of magnitude.

As mentioned, the load F_z = 2.4 is about the minimum load to observe crystallization in the case of frictionless particles.²⁸ Signs of instability are indeed displayed in Fig. 2. The frictionless flowing particles undergo several phase transitions. The system seems steady for a short period (constant K_E in Fig. 2a) with an average volume fraction of about 0.58 (Fig. 2b). Then, the particles rearrange and exceed the random close packing limit, 0.64³⁴ (Fig. 2b). During this first rearrangement, the average kinetic energy per particle slowly increases because the flow becomes asymmetric, with the center of mass of the flowing particles moving towards the bottom plane, as revealed by the fact that its x-velocity becomes progressively more negative (Fig. 2c). At t ≈ 120 000, the frictionless system further re-organizes: the average kinetic energy per particle and the x-velocity of the center of mass reverse their trend, while the average solid volume fraction increases well beyond the random close packing. For a short time interval of almost Δt ∼ 10 000 around t ∼ 140 000, each of the three quantities holds a constant value: K_E ∼ 120, [small nu, Greek, macron] ∼ 0.67 and [small nu, Greek, macron]_x ∼ 0, indicating the existence of a steady, symmetric, crystallized state. This steady state is, however, unstable, and is followed by a further time evolution of the quantities of interest.

Fig. 2 Time evolution of (a) the average kinetic energy per particle, (b) the average volume fraction, (c) the velocity of the center of the mass of the flowing particles in the x-direction, for F_z = 2.4 and different particle friction coefficients.

As the particles become frictional, the crystallization process quickens and a steady state is rapidly attained (Fig. 2a), with the particles moving as a rigid body with either one of the bumpy planes (Fig. 2c), depending on the initial conditions. Only the slightly frictional case (μ_p = 0.01) still shows some sign of instability, with the existence of two possible steady states, one of which is short-lived (60 000 < t < 70 000) at an average solid volume fraction of about 0.64 (Fig. 2b), and one which appears more stable (t > 75 000) at an average solid volume fraction well beyond the random close packing. At larger values of the inter-particle friction, the steady state seems unique, with temporal fluctuations in the average solid volume fraction that diminish as μ_p increases (Fig. 2b).

Fig. 3 shows the profiles of solid volume fraction, flow velocity in the flow direction and granular temperature during the steady state, for F_z = 2.4, when the particles are slightly frictional (μ_p = 0.01) and highly frictional (μ_p = 0.3). The flow height at steady state is H ∼ 12 for both values of μ_p. In the case of slightly frictional particles, we observe a crystallized region in which the solid volume fraction almost reaches the maximum packing for identical spheres, i.e. ν = 0.74 for face-centred cubic (fcc) and hexagonal-closed packed (hcp) structures, in the core region. Shearing is strongly localized in one shear band near the upper plane (z/H > 0.8), whereas the particles in the crystallized region move with an almost uniform velocity (Fig. 3b). The granular temperature is larger in the shear band (Fig. 3c). Similar results were obtained in dense granular shear flows by Alam and Luding³⁵ and Mokshin and Barrat.²³ In the shear band, the solid volume fraction is as low as 0.2 (Fig. 3a). We measure unexpectedly low values of the solid volume fraction near the bottom boundary, where the particles are crystallized, but we deem it an artifact of coarse-graining near a boundary.

Fig. 3 Local time-averaged profiles of (a) solid volume fraction, (b) velocity in the flow direction and (c) granular temperature, for F_z = 2.4 with μ_p = 0.01 (solid lines) and μ_p = 0.3 (dashed lines).

In the case of highly frictional particles, the shear band disappears, and all particles move as a rigid body at the velocity of the lower bumpy plane (Fig. 3b), while experiencing full slip with the upper plane. The granular temperature is at least one order of magnitude less than in the slightly frictional case (friction suppresses the particle agitation) and the solid volume fraction exceeds the random close packing only near the bumpy planes, at z/H ∼ 0.2 and 0.8.

A more detailed picture of the micro-structures is revealed by the time-averaged local order parameters. The pairs for each flowing particle in the system are displayed in Fig. 4a for the slightly frictional case, coloured according to their position in the x−z plane of the domain (Fig. 4b). The crystallized particles are observed to be both in the fcc and hcp phases, but the fcc structures dominate. We notice that the flowing particles near the bottom boundary also align themselves with the boundary particles in a way that yields fcc layering; however, this alignment is not stable over time. The particles in the top shear band region are in between the crystalline state and the fluid state (Fig. 4a).

Fig. 4 (a) Time-averaged local order parameters for flowing particles of μ_p = 0.01 when F_z = 2.4, coloured by their position in the z-direction. The sketched regions for bcc, hcp, fcc, and fluid phases are taken from Eslami et al.³³ (b) A snapshot of the particle positions projected on the x−z plane (particle size is scaled by 50% to facilitate the view).

As the particles become more frictional, the ordering in the core region of the system is destroyed, and most of the particles there belong to the fluid phase (Fig. 5a and b). This disordered region percolates intermittently to the boundaries, although most of the particles adjacent to the bumpy planes remain in the fcc or hcp phase. The crystallized flowing particles align with the boundary particles and form hcp layering, which is stable over time.

Fig. 5 (a) Time-averaged local order parameters for flowing particles of μ_p = 0.3 when F_z = 2.4, coloured by their position in the z-direction. The sketched regions for bcc, hcp, fcc, and fluid phases are taken from Eslami et al.³³ (b) A snapshot of the particle positions projected on the x−z plane (particle size is scaled by 50% to facilitate the view).

Increasing the vertical load by four orders of magnitude results in quicker self-organization of the flowing particles into a crystalline state, as revealed by the average solid volume fraction rapidly increasing to about 0.7, with only a mild dependence on the inter-particle friction coefficient (Fig. 6b). Compared to the small load case, larger fluctuations in the average kinetic energy per particle are observed (Fig. 6a). In most cases, the steady-state configuration is not symmetric with respect to the z-direction, as revealed by the non-vanishing x-velocity of the center of mass of the particles (Fig. 6c).

Fig. 6 Time evolution of (a) the average kinetic energy per particle, (b) the average volume fraction, (c) the velocity of the center of the mass of the flowing particles in the x-direction, for F_z = 24 000 and different particle friction coefficients.

In the case of high load, there is a shear band characterized by low solid volume fraction and high granular temperature localized inside the flow domain (Fig. 7). In the steady state, regardless of the friction coefficient of the flowing particles, the shear band has a thickness of about 2 to 3 diameters. The location of the shear band (i.e., high-temperature region) can be either in the middle of the domain or lean towards one of the bumpy planes (Fig. 7c). The shear band is squeezed between two blocks of crystallized particles moving as rigid bodies with bumpy planes (Fig. 7b).

Fig. 7 Local time-averaged profiles of (a) solid volume fraction, (b) velocity in the flow direction and (c) granular temperature, for F_z = 24 000 with μ_p = 0.01 (solid lines) and μ_p = 0.3 (dashed lines).

The majority of the crystallized particles under high load occupy the fcc phase, whereas a few are observed in the hcp phase, irrespective of the particle friction coefficient (Fig. 8a and 9a for μ_p = 0.01 and 0.3, respectively). The crystallized flowing particles are also perfectly aligned with the boundaries as seen in Fig. 8c and 9c. Even inside the shear bands, the particles are far from the fluid state. The dominance of fcc structures in our results agrees with previous works on the crystallization of colloidal suspensions³⁶ and granular spheres subjected to shearing or vibration,^14,20,37 where it has been suggested that the fcc structure has more stable mechanical properties than the hcp structure. Bai et al.¹⁸ have instead observed that the hcp structure was more favoured than the fcc in their shearing flows, probably because of the curvature of the boundary in their system.

Fig. 8 (a) Time-averaged local order parameters for flowing particles of μ_p = 0.01 when F_z = 24 000, coloured by their position in the z-direction. The sketched regions for bcc, hcp, fcc, and fluid phases are taken from Eslami et al.³³ (b) A snapshot of the particle positions projected on the x−z plane (particle size is scaled by 50% to facilitate the view).

Fig. 9 (a) Time-averaged local order parameters for flowing particles of μ_p = 0.3 when F_z = 24 000, coloured by their position in the z-direction. The sketched regions for bcc, hcp, fcc, and fluid phases are taken from Eslami et al.³³ (b) A snapshot of the particle positions projected on the x−z plane (particle size is scaled by 50% to facilitate the view).

Finally, we plot in Fig. 10 the macroscopic friction, μ, measured in the steady state as a function of the particle friction coefficient, μ_p, in the case of both small and high normal loads. The fact that the crystalline structures are independent of μ_p under high loads and always perfectly align with the bumpy boundaries results in a macroscopic geometrical friction that is also independent of the particle friction and roughly equal to 0.2. When the load is slightly above the minimum to observe crystallization in the case of frictionless particles, we have shown that particle friction controls the transition of a significant portion of the particles into a disordered, fluid state. Unlike solids, the shear resistance of fluids vanishes at zero shear rate. Hence, the unexpected consequence that the macroscopic friction for highly frictional particles is half of that for frictionless particles (Fig. 10) is likely due to the increased portion of the particles occupying the fluid state.

Fig. 10 Macroscopic friction as a function of the inter-particle friction coefficient in the steady state for F_z = 2.4 (circles-solid line) and F_z = 24 000 (squares-dashed line).

4. Conclusions

We have performed discrete element simulations of identical, inelastic spheres sheared, in the absence of gravity, between two parallel, bumpy planes, under a constant load high enough to induce crystallization. The bumpiness of the planes was due to a monolayer of identical, frictionless spheres glued together and arranged in a regular, hexagonal lattice. We have probed the effect of the friction coefficient of the flowing particles on the self-assembly process and the resulting microstructures, for two different loads.

For the smaller load, which is slightly above the minimum to induce crystallization, and frictionless particles, the ordering progresses slowly, and the system does not reach a stable steady state. Switching from frictionless to frictional particles causes quicker self-organization, probably due to the enhanced energy dissipation. For slightly frictional particles, we have observed one shear band adjacent to one of the boundaries and one crystallized region on the other side where all particles move as a plug. The face-centered cubic and hexagonal-closed packed phases coexist in the ordered region, with the fcc dominating, and align with the lattice of the bumpy boundary. Increasing the friction coefficient destroys the ordering in the core region, where many particles experience a phase transition to a fluid phase. The striking result is that the overall shear resistance diminishes by a factor of two as the flowing particles become more frictional.

The particle friction coefficient seems not to affect the flow dynamics and the microstructures as the vertical load increases by four orders of magnitude. The steady state for different particle friction exhibits similar behaviour, with one shear band of thickness 2 to 3 particle diameters in the core region, surrounded by two blocks of particles in the fcc phase moving as rigid bodies with the boundaries. The particles in the ordered regions always form fcc layering with the bumpy boundaries. The overall shear resistance is similar to the case of frictionless particles at smaller load, and roughly independent of the friction coefficient.

It is worthwhile mentioning that the solid volume fraction in our systems is large enough that even in the presence of an interstitial viscous fluid, the inter-particle contact forces would be dominant.³⁸ We therefore expect the present findings to also apply to dense, non-Brownian suspensions.

Our study indicates that the shear resistance of crystallized spheres depends on two competing mechanisms at work: (i) the disordering effect of inter-particle friction, and (ii) the ordering induced by the applied load. If the applied load is just high enough so that order is slightly favoured, then the friction controls the transition of a portion of the particles into a fluid phase, with a significant reduction of the overall shear resistance. This result highlights the potentiality of using granular materials as effective lubricants in view of industrial applications.

Author contributions

All authors designed research. EK performed simulations, evaluated the data and wrote the manuscript. EK and DV post-processed the data. All authors discussed the results, read, revised and approved the final manuscript.

Data availability

Data for this article are available at Zenodo at: https://doi.org/10.5281/zenodo.14609787.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This project has received funding from the HORIZON-EIC-2021-PATHFINDEROPEN-01 N. 101046693, SSLiP project. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or European Commission. Neither the European Union nor the granting authority can be held responsible for them.

Notes and references

1. H. M. Jaeger, S. R. Nagel and R. P. Behringer, Rev. Mod. Phys., 1996, 68, 1259–1273.
2. S. Mandal and D. V. Khakhar, Phys. Rev. E, 2017, 96, 050901.
3. J. C. Tsai, G. A. Voth and J. P. Gollub, Phys. Rev. Lett., 2003, 91, 064301.
4. J. C. Tsai and J. P. Gollub, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 031303.
5. E. R. Nowak, J. B. Knight, M. Povinelli, H. M. Jaeger and S. R. Nagel, Powder Technol., 1997, 94, 79–83.
6. P. Philippe and D. Bideau, Europhys. Lett., 2002, 60, 677.
7. K. E. Daniels and R. P. Behringer, Phys. Rev. Lett., 2005, 94, 168001.
8. M. P. Ciamarra, A. Coniglio, D. D. Martino and M. Nicodemi, Eur. Phys. J. E: Soft Matter Biol. Phys., 2007, 24, 411–415.
9. O. Carvente and J. C. Ruiz-Suárez, Phys. Rev. Lett., 2005, 95, 018001.
10. A. B. Yu, X. Z. An, R. P. Zou, R. Y. Yang and K. Kendall, Phys. Rev. Lett., 2006, 97, 265501.
11. X. Z. An, R. Y. Yang, R. P. Zou and A. B. Yu, Powder Technol., 2008, 188, 102–109.
12. X. Z. An, C. X. Li, R. Y. Yang, R. P. Zou and A. B. Yu, Powder Technol., 2009, 196, 50–55.
13. C. X. Li, X. Z. An, R. Y. Yang, R. P. Zou and A. B. Yu, Powder Technol., 2011, 208, 617–622.
14. D. P. Shinde, A. Mehta and G. C. Barker, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 022204.
15. D. A. Morales-Barrera, G. Rodriguez-Gattorno and O. Carvente-Muñoz, Phys. Rev. Lett., 2018, 121, 074302.
16. R. Amirifar, K. Dong, Q. Zeng, X. An and A. Yu, Powder Technol., 2021, 380, 47–58.
17. D. S. Grebenkov, M. P. Ciamarra, M. Nicodemi and A. Coniglio, Phys. Rev. Lett., 2008, 100, 078001.
18. J. Bai, J. Li, G. Hong, J. Pan and H. Fei, Powder Technol., 2023, 426, 118615.
19. D. Berzi, Phys. Rev. Fluids, 2024, 9, 034304.
20. A. Panaitescu, K. A. Reddy and A. Kudrolli, Phys. Rev. Lett., 2012, 108, 108001.
21. F. Rietz, C. Radin, H. L. Swinney and M. Schröter, Phys. Rev. Lett., 2018, 120, 055701.
22. N. Duff and D. J. Lacks, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 031501.
23. A. V. Mokshin and J.-L. Barrat, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 021505.
24. L. E. Silbert, G. S. Grest, S. J. Plimpton and D. Levine, Phys. Fluids, 2002, 14, 2637–2646.
25. V. Kumaran and S. Maheshwari, Phys. Fluids, 2012, 24, 053302.
26. V. Kumaran and S. Bharathraj, Phys. Fluids, 2013, 25, 070604.
27. G. Yang, S. Zhang and P. Lin, Phys. Fluids, 2022, 34, 013313.
28. D. Vescovi, A. S. de Wijn, G. L. Cross and D. Berzi, Soft Matter, 2024, 20, 8702–8715.
29. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun, 2022, 271, 108171.
30. A. Thornton, T. Weinhart, S. Luding and O. Bokhove, Int. J. Mod. Phys. C, 2012, 23, 1240014.
31. T. Weinhart, A. Thornton, S. Luding and O. Bokhove, Granular Matter, 2012, 14, 289–294.
32. P. Steinhardt, D. Nelson and M. Ronchetti, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 784–805.
33. H. Eslami, P. Sedaghat and F. Müller-Plathe, Phys. Chem. Chem. Phys., 2018, 20, 27059–27068.
34. S. Torquato, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 3170–3182.
35. M. Alam and S. Luding, Phys. Fluids, 2003, 15, 2298–2312.
36. L. V. Woodcock, Nature, 1997, 385, 141–143.
37. S. Heitkam, W. Drenckhan-Andreatta and J. Fröhlich, Phys. Rev. Lett., 2012, 108, 148302.
38. R. Seto, R. Mari, J. F. Morris and M. M. Denn, Phys. Rev. Lett., 2013, 111, 1–5.

---

Footnotes

† https://www.lammps.org

‡ https://www.mercurydpm.org

---