A Computational Study on the Effect of Initial Powder Packing Configuration on Final Sintered Part Microstructures

Sintering, a process used for creating solid components from powders, is widely used in manufacturing. Though sintering is used for metals, polymers, ceramics and composites, it is a preferred manufacturing technique for ceramic components due to their brittleness and high melting temperatures. Sintering is a solid-state process; that is, the powder particles are fused into a bulk material through the application of heat (and generally pressure), without melting. Furthermore, this process densifies the material, removing porosity as the particles fuse together. Thus the final microstructure of the part is highly dependent on the sintering temperature, pressure and duration, as well as the initial packing configuration of the powder particles in the mold.

For ceramics, near-zero porosity is critical for the preservation of desirable mechanical properties: pores reduce the cross-sectional area of a part, increasing the stress concentration it experiences for a given loading state, resulting in nucleation sites for cracks. Refined grain microstructures are also generally desired for high strength, as prescribed by the Hall-Petch effect. Optimizing these properties can be difficult, since the long sintering times needed for full consolidation often produce large grains¹.

Due to financial and time constraints associated with experimental work, it is desirable to design a high-throughput computational model for predicting part characteristics as a function of several input parameters (e.g. powder characteristics and sintering time). This methodology uses existing high-performance computing (HPC) resources to make connections between the initial powder characteristics and the final part microstructure. Accordingly, it is necessary to accurately and efficiently model both the initial packing of the powder in the mold—to form the green body compact—and the subsequent sintering of the green part into a densified component.

Generating synthetic green body compacts remains a computational challenge, especially for non-monomodal particle size distributions. As the size difference between the largest and smallest particles increases, the total number of particles rises sharply and dramatically increases the computational cost of simulating interparticle interactions, even when a fairly simple contact model is used. A packing algorithm recently employed by Srivastava et al.² yielded estimates for packing efficiency in bidisperse (two particle sizes present in the simulation) frictional arrangements of spheres. The maximum packing efficiency achieved in this work was reported to be approximately 0.87 for frictionless spheres with size dispersions up to 40:1, occurring when the cumulative volume fraction (CVF) of small particles in the mixture was about 0.26. Packing efficiencies for spheres with sliding friction (μ_s=0.30) were reduced slightly, showing a peak value of around 0.85. Subsequently, a frictionless model³ was extended to power-law size distributions using an improved neighbor binning algorithm, managing to simulate powder packings with size dispersions up to 100:1.

Starting from the green body compact, the sintering process can be modeled using methods derived from Monte Carlo simulation⁴. These simulations have been parameterized to match experimental findings for copper powder sintering⁵, showing promise for direct application to higher-temperature ceramic sintering processes; furthermore, these Monte Carlo simulations have been used to predict the strain within sintered components that arises from vacancy annihilation⁶.

In this work, both a discrete-element powder packing model and a Monte-Carlo-based sintering model were implemented. In order to generate accurate results without unnecessary computational expense, convergence tests were conducted to determine the appropriate simulation cell size, timestep size, and total number of timesteps for the powder packing model. The discrete-element model was also used to investigate the effects of particle size ratio, cumulative volume fraction and sliding friction coefficient on the packing efficiency in systems with bidisperse and tridisperse size distributions. For the sintering model, the number of timesteps and the discretization mesh density were varied in order to establish a relationship between the two parameters.

Discrete Element Modeling

To understand the impact of powder characteristics (e.g., particle size distribution) on the initial particle packing configuration within a mold, a granular material simulation was constructed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)⁷ developed by Sandia National Laboratories. The simulation employed in this work—a frictional sphere packing model—idealizes all powder particles as elastic spheres which experience friction against the other particles and the walls of the container (wherever non-periodic boundaries exist).

Elastic interactions between particles were simulated according to a Hertz-Mindlin contact model^8,9 via the LAMMPS hertz/material command. This model computes the force between particles i and j as

where k_n is the elastic stiffness of the material, R_eff=R_i R_j⁄(R_i+R_j) is the effective interaction radius, and is the particle overlap. Here, R_i and R_j are the radii of particles i and j, and is the separation between the particle centers. The direction vector for the interparticle force is

In addition to elastic interactions, the Hertz-Mindlin contact model includes particle-particle and particle-container frictional interactions. For this work, only a nonzero tangential sliding friction was simulated, with a friction coefficient μ_s. Only non-periodic simulation cell boundaries were modeled with frictional interactions; for periodic boundaries, no particle-container interaction was present.

Initially, the powder packing configurations were determined via a granular pouring model. This approach initiates particles at the top of a vertically elongated simulation cell and allows them to fall to the bottom due to gravity. As more particles fall, they pile up at the bottom of the cell and eventually settle into a stable configuration. Once a sufficient number of particles has been

introduced, the simulation volume is truncated and a small number of additional timesteps is allowed to elapse. Then the lower 75% of the simulation volume is retained in order to isolate the most densely packed portion of the overall configuration, and this reduced volume is taken as the final powder packing configuration.

Figure 1. Visualization of the granular pouring model used in this work.

Following the work of Srivastava et al.² and Santos et al.¹⁰ a confinement-packing model was also implemented. This approach uses a very large initial simulation volume, populated randomly with particles, wherein overlapping particles are deleted and a small confinement pressure is applied to jam the particles together.

A 30R×30R×30R initial simulation volume was selected in order to produce a roughly 10R×10R×10R final volume, where R is the radius of the largest powder particle in the mixture. The applied confinement pressure was chosen to be

where R_min is the radius of the smallest particle in the simulation. For a mono-sized set of particles, R_min=R. It should be noted that the state of applied stress is purely hydrostatic; that is, P_a,xx=P_a,yy=P_a,zz=P_a and P_a,xy=P_a,yx=P_a,yz=P_a,zy=P_a,zx=P_a,xz=0.

Figure 2. Visualization of the confinement-packing model used in this work.

Determination of Appropriate Simulation Cell Size

Due to the use of periodic boundary conditions in the LAMMPS simulation, an insufficiently large simulation volume can produce direct or indirect particle self-interaction, rendering the simulation unstable. To mitigate this effect, a convergence study was performed to determine the appropriate simulation cell size. This study entailed simulations of the same powder pouring process, but with successively larger simulation cells in the x- and y-dimensions. The end-of-run packing efficiency, ϕ, was computed as an average of 10 simulations for each cell size; the mean final packing efficiency, μ_ϕ, and standard deviation in the final packing efficiency, σ_ϕ, were used to determine convergence. This process was subsequently repeated for the packing-confinement model.

Determination of Appropriate Timestep Size

Following the work of Monti et al.³, the appropriate timestep size was assumed to depend on the smallest particle present in the mixture. If the time resolution of the simulation is too coarse, excessive movement of particle centers over a timestep can cause large elastic deformations of small particles, generating unrealistically high interparticle forces and very high particle velocities. This can prevent the simulation from ever converging to a stable packing arrangement, resulting in anomalously low packing efficiencies. Accordingly, the timestep was chosen to be Δt=ατ, where α is a scaling coefficient (to be determined experimentally) and τ is defined as

where is the effective mass of the smallest particle present in the system, ρ being the mass density of the particles. In order to estimate an appropriate α value for this work, several simulations with different α values were performed; the average and standard deviation of the packing efficiency were examined at each α value to determine adequate convergence. These simulations were performed using the confinement-packing model, for both a monosize distribution of particles and a bidisperse distribution with a 10:1 size ratio.

Determination of Appropriate Number of Timesteps

The acceptable number of timesteps, N_max, was determined using the same approach: a large number of otherwise-identical simulations with different N_max values were carried out, and convergence was gauged using μ_ϕ and σ_ϕ for ten simulations at each N_max value.

Characterization of Powder Packing Configurations

Final powder packing configurations were characterized according to three separate metrics: packing efficiency, lacunarity and Voronoi entropy. The packing efficiency was straightforwardly taken to be the cumulative volume of all particles present at the end of the simulation, normalized by the final volume of the simulation cell:

Here, V_i and R_i are the volume and radius of the i^th particle. Lacunarity, which can be thought of as a measure of inhomogeneity or texture, can be quantified by partitioning the simulation cell into sub-volumes and counting the number of particle centers present within each sub-volume. Using the work of Chekuryaev et al.¹¹ as a guide, the lacunarity parameter Λ was computed according to

where σ is the standard deviation in the number of particle centers found in each sub-volume, and μ is the mean number of particle centers per sub-volume. The Voronoi entropy was also computed for final powder packing configurations, via the following equation:

where P_i is the probability of a given Voronoi polyhedron having i faces. Generation of Voronoi polyhedra¹² was performed automatically in LAMMPS.

Figure 3. Visualization of Voronoi polyhedra within a random sphere-packing arrangement.

Monte Carlo Sintering Model

The sintering process was modeled using a Monte-Carlo-based approach via the Stochastic Parallel PARticle Kinetic Simulator (SPPARKS) simulation package. The particular model developed for this work uses a rejection kinetic Monte Carlo (rKMC) algorithm to model three types of events—grain growth, pore migration and pore annihilation—as the initial ‘green compact’ is sintered into a bulk material.

The green compact modeled at the start of the SPPARKS sintering process is derived from the final powder packing configuration of a LAMMPS simulation. Each spherical particle is assumed to be a single crystal, and is accordingly assigned a random ‘spin’ to represent its crystallographic orientation. The SPPARKS sintering application discretizes the LAMMPS output configuration into a voxelized domain, where voxels are either filled (order parameter larger than 0) or unfilled (0) with material. The density of the discretized domain (or, alternatively, the number of voxels in the simulation cell) therefore has a pronounced effect on the behavior of the model. To characterize this effect, a series of SPPARKS simulations were performed with the spacing between mesh points, Δx, being varied. Each datum shown for these convergence tests is averaged over five SPPARKS runs with identical inputs.

Figure 4. Evolution of part microstructure within SPPARKS sintering simulation. Note the initial increase in exterior roughness due to pore annihilation. Once the porosity has sufficiently decreased, grain growth becomes the dominant process.

In the sintering model, the total system energy is computed according to the work of Tikare et al.⁵:

where N is the total number of sites (or voxels) in the simulation, n is the number of neighbors for each site, δ(q_i,q_j ) is the Kronecker delta applied to the state q_i of each site i and the states q_j of each of its neighbor sites j.

The frequency at which some event i is attempted is given by

where E_i is the activation energy for the event, ω_0,i is a constant prefactor, T is the temperature, and k_B is Boltzmann’s constant. Each attempted process is then either accepted or rejected according to an acceptance probability given by

where ΔE is the global energy change resulting from the proposed event. The probability of event i occurring, then, is given by

Accordingly, the SPPARKS sintering model takes as inputs the attempt rate prefactors for grain growth, pore migration and pore annihilation (ω_0,gg, ω_0,pm and ω_0,pa, respectively) and the event temperatures . For the purposes of the tests performed in this work, the E_i were set to zero, meaning that ω_i=ω_0,i; in the future, these values can be replaced with estimates for a particular material to give reasonable variation in event attempt rates with temperature.

Figure 5. Representation of pore migration event, showing a vacancy break away from a pore by swapping positions with a filled mesh site. Different-colored regions represent distinct grains.

Figure 6. Representation of vacancy annihilation event, showing a column of sites move through the grain center of mass to fill in the vacancy at the grain boundary.

Figure 7. Representation of grain growth, showing a site at the boundary changing its ‘spin’ to grow the leftward grain at the expense of the rightward one.

In the sintering model, grain growth is proposed as a random alteration of the ‘spin’ of a given filled site to that of one of the neighboring filled sites; the probability of acceptance of a grain growth event is computed as described above. Thus, over a large number of timesteps, the system will evolve to favor certain grains over others, allowing them to grow at the expense of their neighbors. Pore migration in the sintering model is also straightforward: an unfilled site is selected, and if it has a filled neighboring site, then the two sites may be exchanged. The newly filled site has its spin set to the most energetically favored grain orientation. Through this process, unfilled sites are allowed to break off from pores (agglomerations of unfilled sites) to form vacancies. The third process modeled, vacancy annihilation, deals specifically with these isolated unfilled sites. To annihilate a vacancy, a path is defined from the vacancy to the outside surface of the part, passing through the center of mass of the nearest grain. All sites along the path, filled or unfilled, are then moved one position toward the vacancy, annihilating it and creating a vacant site at the surface of the part.

LAMMPS Sphere Packing Model Results

While close packing of uniformly-sized spheres produces a theoretical packing efficiency of approximately 0.74, random (frictionless) packing is understood to produce a packing efficiency of only around 0.64¹². For the granular pouring model, the converged value was seen to be approximately 0.624, which is lower than the theoretical maximum packing efficiency due to the inclusion of frictional interactions between particles. The end-of-run packing efficiency was found to deviate wildly from this anticipated value when the minimum simulation cell dimension was below about 6R; beyond 10R, the results were deemed to be sufficiently converged, as the standard deviation in packing efficiency between runs dropped to below 1%. Additionally, the variation between runs—arising from the random nature of the simulation—accounted for more variability than did the cell size, so long as the simulation cell was no smaller than about 10R in any dimension.

Figure 8. Convergence of packing efficiency with respect to simulation cell size in the granular pouring model. Each value is averaged from 10 independent simulations.

Figure 9. Standard deviation in packing efficiency over varying simulation cell sizes for the granular pouring model. Each datum is calculated from 10 independent simulations.

 

Figure 10. Variation in packing efficiency for the confinement-packing model with a bidisperse particle size distribution and a 5:1 size ratio.

Figure 11. Standard deviation in packing efficiency at different simulation box sizes for the confinement-packing model. Statistics for each datum are generated from the outputs of 10 simulations.

For the timestep convergence, the acceptable timestep multiplier value was found to be around α=0.10; beyond this value, the mean packing efficiency settles to a value between 0.6291 and

0.6302, with the standard deviation in final packing efficiency between identical runs falling to about 0.002.

Figure 12. Convergence of packing efficiency with respect to timestep size in the confinement-packing model. The model is seen to give stable results even with large timesteps.

Figure 13. Variation in packing efficiency for a bidisperse particle size distribution with a 10:1 size ratio. Each datum represents a single simulation result.

For comparison, both the pouring model and the confinement-packing model were examined for convergence over a wide range of N_max values. Using a very conservative timestep (Δt=0.01τ), it is evident that the packing-confinement model converges very quickly to a stable packing efficiency: after approximately one million timesteps, or 1×10⁴ τ, ϕ has converged to approximately 0.624 (see Figure 14). By contrast, the pouring model fails to converge after 10^8 timesteps, or 1×10^6 τ. Accordingly, the packing-confinement model was deemed superior to the pouring model, and the remainder of the LAMMPS simulation work was performed using the more efficient methodology.

Figure 14. Evolution of final packing efficiency with respect to total runtime of the granular pouring model (top), and the confinement-packing model (bottom). The confinement-packing model can be seen to exhibit much faster convergence, even with a very conservative timestep size (Δt=0.01τ). Both tests shown are for bidisperse frictional packings with a 2:1 particle size ratio.

Figure 15. Variation in packing efficiency and Voronoi entropy with respect to the number of allowed simulation timesteps. Note that greater particle size dispersion contributes to both increased packing efficiency and lower (smaller magnitude) Voronoi entropy.

As indicated by prior work on bidisperse packing of spheres2, the packing efficiency increases substantially as the size dispersion increases. Furthermore, an increase in the friction coefficient naturally reduces packing efficiencies. A corresponding reduction in the magnitude of the Voronoi entropy is seen as the simulations converge.

Even with the timestep scaled according to the effective mass of the smallest particle in the system (see Materials & Methods section for details), it is clear that a larger size dispersion contributes significantly to the computational cost of the simulation. The packing efficiency data in Figure 15 show the 5:1 size ratio taking roughly tenfold more timesteps than the 2:1 size ratio to converge, representing an enormous increase in resource consumption. Hence the ability of the confinement-packing model to converge rapidly even with a conservative timestep yields a significant advantage over the pouring model.

Figure 16. Comparison of packing efficiencies for various particle size dispersions, cumulative volume fractions, and sliding friction coefficients. Lower friction coefficients and higher size dispersions result in increased packing efficiency. The optimal CVF of small particles is always roughly 0.3, although it shows some variation.

Figure 17. Comparison of packing efficiency results from Srivastava et al.2 and the current work for bidisperse distributions of frictional elastic spheres with μ_s=0.30.

The lacunarity parameter was also computed for many of these packed configurations, showing the inhomogeneity of the particle mixture. As is evident from Figure 18, the computed lacunarity is strongly dependent on the size of the sub-volumes into which the sample is divided. When the sub-volume is a very small portion of the simulation cell, the standard deviation in the number of particle centers per sub-volume is large compared to the mean value, so the lacunarity parameter increases in magnitude rapidly. At the other extreme, if the simulation cell contains only one sub-volume, then the lacunarity is zero. This size variation also means that there are variable measures of lacunarity that depend on which particles are included in the analysis. Figure 18 shows lacunarity values for the same structure—a 5:1 bidisperse arrangement of randomly packed spheres—computed with different sub-volume sizes and counting choices. If only the large particles are counted, then the lacunarities are quite high, indicating a rough texture. However, if only the small particles are counted, far greater homogeneity is seen.

Figure 18. Variation in lacunarity parameter Λ with sub-volume edge length, x (as a multiple of the large particle radius R). The data is presented as a linear-linear plot (top) and as a log-log plot (bottom), showing the different values obtained from different particle center counting methods.

In addition to various bidisperse particle size distributions, a tridisperse size distribution with a 4:2:1 size ratio was modeled. As can be seen in Figure 19, the sides of the ternary diagram are consistent with the bidisperse packing efficiency results. In the case of this size distribution, the maximum packing efficiency lies along the rightward side of the ternary (which represents an 8:1 bidisperse mixture). It is possible that, with greater size ratios, the optimal packing efficiency may shift away from the sides of the ternary; however, this was not pursued in the current work.

Figure 19. Packing efficiency of powder with a tridisperse size distribution and a 4:2:1 size ratio. The terminal labels S, M and L indicate the smallest, middle and largest monosize particle data, respectively.

The convergence behavior of the SPPARKS model was seen to exhibit a strong dependence on the number of voxels in the simulation cell. Figure 20 shows the density of the part (relative to a bulk material) increasing over the course of the sintering simulations with different Δx values. It is evident from these results that an increase in the resolution of the discretization mesh dramatically increases the cost of the simulation. This is intuitive, as a halving of the mesh spacing results in an eightfold increase in the total number of mesh sites. If the mesh is too coarse (corresponding to a spacing of around Δx≥0.25), then significant variation appears between individual runs and the results become less stable.

Figure 20. SPPARKS sintering curves for different discretization mesh spacings, Δx. Note that the denser meshes, having considerably more voxels, require a far larger number of timesteps to converge. If the discretization is too coarse, the results become unstable.

The convergence speed with respect to the number of simulation timesteps was initially assumed to trend roughly linearly with the total number of voxels, meaning that the mesh spacing Δx would have a cubic relationship with the required number of steps. To examine this possibility, the sintering curves were plotted against the number of timesteps, N, scaled by the voxel volume, Δx^3. From the data in Figure 21, it is clear that this relationship is not precisely linear, since the sintering curves do not appear to converge to a stable value as Δx is reduced. This implies that, while the number of required steps is somewhat volume-dependent, there are other dependencies that alter the overall behavior of the model. This is perhaps unsurprising, as processes such as grain growth will exhibit area dependence.

Figure 21. Sintering curves scaled by voxel volume, showing similar behavior but convergence to slightly different final relative densities.

As is evident in Figure 22, the behavior of the model becomes more predictable as the mesh spacing is decreased. However, the end-of-run relative density of the sintered part appears to converge to slightly different values as the mesh size is altered. So, although not fully converged to a predictable final density, the randomness of the model is greatly reduced by refinement of the discretization mesh.

Figure 22. Standard deviation in final relative densities for different mesh spacings, showing improvement in simulation stability with mesh refinement.

The fact that the final relative density, ρ_f^rel, does not smoothly converge to a stable value as the mesh size is decreased indicates different dependences of the kinetic processes modeled. In order to precisely determine the behavior of the model overall, the volume and area dependencies of the different kinetic processes will need to be independently quantified.

Both of the LAMMPS powder packing models explored in this work (granular pouring and confinement-packing) were seen to produce realistic final packing configurations close to theoretical packing efficiencies for monosized sphere arrangements. The confinement-packing model, however, was found to converge to reasonable estimates of packing efficiency in a far smaller number of timesteps than the granular pouring model, owing to the omission of the actual pouring and settling processes required for the latter model. Thus the confinement-packing model can be used to generate reasonable estimates for packing efficiency, lacunarity and Voronoi entropy without unnecessary resource use, allowing for far greater simulation throughput.

The packing-confinement model was used to generate estimates of powder packing efficiency over a range of friction coefficients, bidisperse particle size ratios, and cumulative volume fractions, and showed reasonably good agreement with prior work. Furthermore, the convergence tests demonstrated the necessity of allowing larger numbers of timesteps in high-size-ratio simulations, accentuating the advantage of the confinement-packing model.

The discretization mesh density in SPPARKS was found to have a dramatic influence on both the stability of the model’s results and the required runtime for the simulation. The total number of mesh sites (or voxels) in the simulation space was found to have a roughly linear effect on the required number of timesteps, although it remains clear that other effects are present.

In order to fully implement the sintering model, it will be necessary to fully characterize the model’s dependence on the relative event probabilities by examining each type of event independently and in greater detail—the grain growth, pore migration and vacancy annihilation processes have different responses to mesh refinement. Furthermore, available experimental data for common materials (e.g. alumina) should be used to give the event attempt rates reasonable temperature dependence; ab-initio or atomistic results may stand in where experimental data is lacking. Completion of these additional tasks will allow for rapid, reliable powder packing and sintering simulations, yielding a high-throughput computational asset for guiding future experiments.

This work was supported in part by High-Performance Computing (HPC) resources from the U.S. Department of Defense (DoD) High-Performance Computing Modernization Program (HPCMP), in collaboration with an appointment to the HPC Internship Program (HIP) administered by the Oak Ridge Institute for Science and Education (ORISE) through an interagency agreement between the U.S. Department of Energy (DoE) and the DoD. ORISE is managed by Oak Ridge Associated Universities (ORAU) under DoE contract number DE-SC0014664. All opinions expressed in this work are the author’s and do not necessarily reflect the policies and views of DOD, DOE, or ORAU/ORISE.

1. Yu Pan, He Li, Yongsheng Liu, Yansong Liu, Kehui Hu, Ning Wang, Zhigang Lu, Jingjing Liang, and Sujun He. 2020. Effect of Holding Time During Sintering on Microstructure and Properties of 3D Printed Alumina Ceramics. Front. Mater. 7, (April 2020), 54. DOI:https://doi.org/10.3389/fmats.2020.00054

2. Ishan Srivastava, Scott A. Roberts, Joel T. Clemmer, Leonardo E. Silbert, Jeremy B. Lechman, and Gary S. Grest. 2021. Jamming of bidisperse frictional spheres. Phys. Rev. Research 3, 3 (August 2021), L032042. DOI:https://doi.org/10.1103/PhysRevResearch.3.L032042

3. Joseph M. Monti, Joel T. Clemmer, Ishan Srivastava, Leonardo E. Silbert, Gary S. Grest, and Jeremy B. Lechman. 2022. Large-scale frictionless jamming with power-law particle size distributions. Phys. Rev. E 106, 3 (September 2022), 034901. DOI:https://doi.org/10.1103/PhysRevE.106.034901

4. Veena Tikare, Michael Braginsky, and Eugene A. Olevsky. 2003. Numerical Simulation of Solid-State Sintering: I, Sintering of Three Particles. Journal of the American Ceramic Society 86, 1 (January 2003), 49–53. DOI:https://doi.org/10.1111/j.1151-2916.2003.tb03276.x

5. Veena Tikare, Michael Braginsky, Didier Bouvard, and Alexander Vagnon. 2010. Numerical simulation of microstructural evolution during sintering at the mesoscale in a 3D powder compact. Computational Materials Science 48, 2 (April 2010), 317–325. DOI:https://doi.org/10.1016/j.commatsci.2010.01.013

6. R. Bjørk, H.L. Frandsen, V. Tikare, E. Olevsky, and N. Pryds. 2014. Strain in the mesoscale kinetic Monte Carlo model for sintering. Computational Materials Science 82, (February 2014), 293–297. DOI:https://doi.org/10.1016/j.commatsci.2013.09.051

7. Steve Plimpton. 1995. Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics 117, 1 (March 1995), 1–19. DOI:https://doi.org/10.1006/jcph.1995.1039

8. R. D. Mindlin. 1949. Compliance of Elastic Bodies in Contact. Journal of Applied Mechanics 16, 3 (September 1949), 259–268. DOI:https://doi.org/10.1115/1.4009973

9. R. D. Mindlin and H. Deresiewicz. 1953. Elastic Spheres in Contact Under Varying Oblique Forces. Journal of Applied Mechanics 20, 3 (September 1953), 327–344. DOI:https://doi.org/10.1115/1.4010702

10. A. P. Santos, Dan S. Bolintineanu, Gary S. Grest, Jeremy B. Lechman, Steven J. Plimpton, Ishan Srivastava, and Leonardo E. Silbert. 2020. Granular packings with sliding, rolling, and twisting friction. Phys. Rev. E 102, 3 (September 2020), 032903. DOI:https://doi.org/10.1103/PhysRevE.102.032903

11. Andrey G. Chekuryaev, Maxim M. Sychov, Sergey N. Perevislov, and Vladimir N. Ulanov. 2023. Digital Characteristics of Microstructure of Diamond—Silicon Carbide Composites. Ceramics 6, 2 (May 2023), 1067–1077. DOI:https://doi.org/10.3390/ceramics6020063

12. Witold Brostow, Jean-Pierre Dussault, and Bennett L Fox. 1978. Construction of Voronoi polyhedra. Journal of Computational Physics 29, 1 (October 1978), 81–92. DOI:https://doi.org/10.1016/0021-9991(78)90110-9

13. Corey S. O’Hern, Leonardo E. Silbert, Andrea J. Liu, and Sidney R. Nagel. 2003. Jamming at zero temperature and zero applied stress: The epitome of disorder. Phys. Rev. E 68, 1 (July 2003), 011306. DOI:https://doi.org/10.1103/PhysRevE.68.011306