Geometrical properties of rigid frictionless granular packings as a function of particle size and shape

TL;DR

This study uses 3D simulations to explore how particle size and shape affect the packing density and microstructure of frictionless granular assemblies, revealing shape and polydispersity influences on solid fraction and order.

Contribution

It provides new insights into the effects of particle shape and polydispersity on granular packing properties, confirming the isostatic conjecture and the impact of particle reorientation.

Findings

Disordered packings of pinacoids achieve higher solid fractions than spherical packings.

Solid fraction increases with bidispersity, peaking at 30% small particles.

Particle shape has a minor effect compared to polydispersity on packing density.

Abstract

Three-dimensional discrete numerical simulation is used to investigate the properties of close-packed frictionless granular assemblies as a function of particle polydispersity and shape. Unlike some experimental results, simulations show that disordered packings of pinacoids (eight-face convex polyhedron) achieve higher solid fraction values than amorphous packings of spherical or rounded particles, thus fulfilling the analogue of Ulam's conjecture stated by Jiao and co-workers for random packings [Y. Jiao and S. Torquato, Phys. Rev. E $\textbf{84}$, $041309$ ($2011$)]. This seeming discrepancy between experimental and numerical results is believed to lie with difficulties in overcoming interparticle friction through experimental densification processes. Moreover, solid fraction is shown to increase further with bidispersity and peak when the volume proportion of small particles reaches $30\%$. Contrarywise, substituting up to $50\%$ of flat pinacoids for isometric ones yields solid fraction decrease, especially when flat particles are also elongated. Nevertheless, particle shape seems to play a minor role on packing solid fraction compared to polydispersity. Additional investigations focused on the packing microstructure confirm that pinacoid packings fulfill the isostatic conjecture and that they are free of order except beyond $30$ to $50\%$ of flat or flat \& elongated polyhedra in the packing. This order increase progressively takes the form of a nematic phase caused by the reorientation of flat or flat \& elongated particles to minimize the packing potential energy. Simultaneously, this reorientation seems to increase the solid fraction value slightly above the maximum achieved by monodisperse isometric pinacoids, as well as the coordination number. Finally, partial substitution of elongated pinacoids for isometric ones has limited effect on packing solid fraction or order.

Tables

Table 1. Table 1: Size (diameter d 𝑑 d for spheres, d = L 2 + G 2 𝑑 superscript 𝐿 2 superscript 𝐺 2 d=\sqrt{L^{2}+G^{2}} for pinacoids) and aspect ratio of simulated particles.

Particles	Size	$L/G$	$G/E$
spheres (large, small)	$d$, $d/3$	$-$	$-$
isometric pinacoids (large, small)	$d$, $d/3$	$1$	$1$
elongated pinacoids	$d$	$2$	$1$
flat pinacoids	$d$	$1$	$3$
flat & elongated pinacoids	$d$	$2$	$3$

Table 2. Table 2: Mean proportion by volume of small X d / 3 subscript 𝑋 𝑑 3 X_{d/3} , elongated X P subscript 𝑋 𝑃 X_{P} , flat X O subscript 𝑋 𝑂 X_{O} , or flat & elongated X P ​ O subscript 𝑋 𝑃 𝑂 X_{PO} particles substituted for large isometric ones in each binary mix simulated (e.g. X P = 14 % subscript 𝑋 𝑃 percent 14 X_{P}=14\% refers to the bidisperse pinacoid packing comprising 14 % percent 14 14\% of elongated pinacoids and 86 % percent 86 86\% of large isometric ones).

Spheres (%)	Pinacoids (%)
$X_{d/3}$	$X_{d/3}$	$X_P$	$X_O$	$X_{PO}$
$0$	$0$	$14$	$12$	$14$
$3$	$13$	$20$	$28$	$29$
$13$	$30$	$36$	$47$	$50$
$23$	$51$	$49$	$77$	$70$
$30$	$64$	$70$	$92$	$100$
$41$	-	$100$	$100$	-
$55$	-	-	-	-

Table 3. Table 3: Summary of simulation parameters. 3000 3000 3000 particles, horizontal dimension of square simulation cell L 𝐿 L normalized by particle size d 𝑑 d (periodic boundary conditions apply along x 𝑥 x and y 𝑦 y directions), friction μ 𝜇 \mu , restitution coefficients e N , T subscript 𝑒 𝑁 𝑇 e_{N,T} and time steps Δ ​ T Δ 𝑇 \Delta T normalized by d / g 𝑑 𝑔 \sqrt{d/g} .

Particle	$L/d$	$\mu$	$e_{N,T}$	$\Delta T/\sqrt{d/g}$
Spheres and	$3$ to $8$	$0$	$0$	$5\times {10}^{-3}$
isometric pinacoids
elongated pinacoids	$8$ to $11$	$0$	$0$	$5\times {10}^{-3}$
flat pinacoids	$6$ to $8$	$0$	$0$	$5\times {10}^{-3}$
flat & elongated	$8$	$0$	$0$	$5\times {10}^{-3}$
pinacoids

Table 4. Table 4: Minimum solid fraction values obtained for packings made of 100 % percent 100 100\% isometric, 100 % percent 100 100\% elongated, 100 % percent 100 100\% flat or 100 % percent 100 100\% flat & elongated pinacoids.

Pinacoid shape	$100\%$	$100\%$	$100\%$	$100\%$
	isometric	elongated	flat	flat & elongated
${\varphi }_{min}$	$0.664$	$0.662$	$0.648$	$0.654$

Equations

$$V=\frac{EG}{2}(L-\frac{E}{3\sqrt{3}}).$$

$$\sum F<10^{-4}{d}^{2}P,$$

$$\sum M<10^{-4}{d}^{3}P,$$

$$E_c<10^{-8}{d}^{3}P,$$

$$\sum _{particles}\rho Vg=F_S.$$

$$k\sum _{particles}\rho Vgz=\sum _{contacts}{F}_C^z{\delta }_C^z$$

$$N(Q_{00}^2)=N(0)+[N(1)-N(0)].[1-\exp (-\frac{Q_{00}^2}{Q_{00,c}^2})]$$

Full text

I Introduction

Granular materials cover a large variety of natural and industrial matter as diverse as soil, rock, building materials, cereals, or drug capsules. Controling their density, defined as the volume fraction occupied by grains (also called solid fraction), is critical to study their mechanical behavior - which may be gaslike, liquidlike or solidlike - and reduce their manufacturing, storage, or packaging costs. Over the centuries, scientists and engineers have attempted to predict the maximum solid fraction achievable by solidlike assemblies - * or packings* - of hard particles knowning their individual geometrical characteristics.

Due to their apparent simplicity, monodisperse hard sphere packings have received wide interest and various outstanding solid fraction values have been reported, ranging from $\phi=0.555\pm 0.005$ for random loose packings Onoda and Liniger (1990) up to $\phi=\pi/\sqrt{18}\approx 0.74$ for crystal-ordered packings Hales (2005). In the former case, the packing structure is often stabilised by interparticle friction in a hypostatic state Makse et al. (2000), that is with less particle constraints than their number of degrees of freedom, whereas the latter case seems to be favored by long term cyclic shear Panaitescu et al. (2012) and refers to highly hyperstatic structures. In between these two extremes lies the random close packing (RCP) state, equivalently defined as the maximally randomly jammed state Torquato et al. (2000) or stable equilibrium state under isotropic pressure of frictionless particles devoid of crystal nucleus Roux (2004). Such a state may be repeatedly achieved following experimental protocols, e.g. gently kneading waxed balls enclosed in a rubber membrane Bernal and Mason (1960), pouring particles at a controlled rate into a cylinder Macrae and Gray (1961), or vertically shaking them within a container Scott and Kilgour (1969), as well as numerically using purely geometrical Lubachevsky and Stillinger (1990) or mechanically-based protocols Makse et al. (2000); Silbert et al. (2002). Whatever the protocol, the generally agreed solid fraction of sphere packings in RCP state is $\phi_{RCP}=0.6366\pm 0.0005$ Scott and Kilgour (1969); Cumberland and Crawford (1987); O’Hern et al. (2002); Silbert et al. (2002) and their structure is known to be isostatic with a mean number of contacts per spheres—or $coordination$ $number$—equal to $z=6$.

However, real granular media are very seldom made of spheres. In particular, natural $aggregates$—i.e., granular materials for construction and civil engineering uses—may at best be sphere-like (rounded) when extracted from alluvial deposit, but most of them are more or less convex polyhedron-like (angular) as a consequence of their processing from the crushing of massive rock deposits. Unfortunately, literature on assessing the packing solid fraction of convex nonspherical particles is far less abundant. Yet such particles would have higher maximum solid fraction values than spheres according to Ulam’s conjecture Gardner (2001). Indeed, Ulam’s conjecture was verified numerically for dense crystal packings of spheroids of axes ratio larger than $\sqrt{3}$ Donev et al. (2004a), regular tetrahedra dimers Chen et al. (2010), Platonic and Archimedean solids Betke and Henk (2000); de Graaf et al. (2011), whose reported maximum solid fraction values are respectively $\phi_{spheroids}\approx 0.7707$, $\phi_{tetra}\approx 0.8663$ and $\phi_{Pl\&Ar}$ in the range $0.784$ for truncated icosahedra to $1$ for cubes. Ulam’s conjecture was also investigated experimentally for rounded frictionless particles such as ellipsoids of axes ratio $1.25:1:0.8$ and angular grains such as tetrahedral dice. Surprisingly, both particle shapes were found able to pack randomly as densely as crystal-ordered spheres, with solid fraction values of $\phi=0.74\pm 0.005$ for the former Man et al. (2005), and an even denser $\phi=0.76\pm 0.02$ for the latter Jaoshvili et al. (2010) in the limit of infinite system size. In fact, additional simulations revealed that dense amorphous packings of frictionless non-cubical Platonic or axisymmetric-low-aspect-ratio particles would achieve solid fraction values ranging between those of spheres in the RCP and the dense crystal states Donev et al. (2004b); Baule et al. (2013); Jiao and Torquato (2011). This conclusion led Jiao and co-workers to suggest that, among convex particles with moderate asphericity (low aspect ratio), spheres possess the lowest solid fraction in RCP state, which is known as the analogue of Ulam’s conjecture for random packings Jiao and Torquato (2011). According to microstructural investigations of numerical packings, the high level of solid fraction observed with nonspherical particles is caused by the higher coordination number needed to constrain their additional rotational degrees of freedom and achieve jamming Donev et al. (2004b). Note however that amorphous packings of frictionless Platonic solids are isostatic Jaoshvili et al. (2010); Jiao and Torquato (2011), whereas axisymmetric particles were found hypostatic at least for aspect ratio smaller than $1.5$ Donev et al. (2004b, 2007); Baule et al. (2013).

In order to increase further the representativeness of real granular media in packing solid fraction studies, one shall account for the particle size distribution which is seldom monodisperse, and consider higher particle aspect ratios. Indeed, from field experience with aggregates, polydispersity increases the packing solid fraction Caquot (1937); Powers (1968) and bidisperse rounded particles tend to pack denser than angular ones Sedran et al. (1994). More or less empirical models have been designed to mimic these phenomena Stovall et al. (1986); Cumberland and Crawford (1987); Ouchiyama and Tanaka (1989); de Larrard (1999); Liu and Ha (2002). A few numerical studies have been published Roux et al. (2007); Farr and Groot (2010), which confirm higher solid fraction values of bidisperse sphere packings compared to monodisperse ones (up to $0.827$ for a particle diameter ratio of $10$), but little to no similar study dealing with polydisperse polyhedron packing was identified. More literature examines the role played by large particle aspect ratio on packing density. Indeed, frictionless axisymmetric particles with length to diameter ratio larger than $2$ to $3$ were reported to exhibit smaller RCP solid fractions than spheres Kyrylyuk and Philipse (2011); Baule et al. (2013) as a consequence of increased excluded volume effects.

For both sphere and convex nonspherical particle packings, numerical approaches proved essential to investigate their microstructure and confirm their amorphous nature. Indeed, the presence of translational order may easily be checked using common statistical tools such as the $pair$ $correlation$ $function$ (or $radial$ $distribution$ $function$), which describes the probability of finding a pair of particles a distance $r$ apart relative to the probability expected for a completely random distribution at the same density Allen and Tildesley (1989). According to this definition, the pair correlation function of a perfectly amorphous packing is $1$ regardless of $r$ value. The pair correlation function of $dense$ amorphous sphere packings obtained from numerical simulations classically tends towards $1$ beyond a few particle diameters Silbert et al. (2002); Agnolin and Roux (2007), thus evidencing the absence of long range translational order. Similar results evidencing even less translational order were reported from nonspherical particle packings Williams and Philipse (2003); Jaoshvili et al. (2010), as a consequence of their loss of rotational symmetry. For the latter, further checking is needed to ensure the absence of orientational order, which may be achieved using the $nematic$ $order$ $parameter$ and the $biaxial$ $parameter$ John et al. (2004); Camp and Allen (1997). These parameters assess the level of alignment of respectively one and two of the particle inertia axes, and they vary between [math] (no alignment) and $1$ (perfect alignment). Several authors scrutinizing the onset of orientational order in initially random convex nonspherical particle packings subjected to a geometrical densification process have reported amorphous-nematic transition characterized by a nematic order parameter jump in the range $[0.2;0.5]$ Allen et al. (1993); Camp and Allen (1997), whereas the nematic-biaxial transition occured for values of the biaxial parameter in excess of $0.2$ Camp and Allen (1997).

The present paper explores the influence of size, angularity and aspect ratio of frictionless particles on their ability to achieve maximally randomly jammed packings. Further to determining maximum solid fraction values, our goal is to understand how the packing microstructure is affected by particle geometrical characteristics. For this purpose, numerical simulations have been performed so as to mimic a mechanically-based densification process, thus allowing comparisons with similar laboratory experiences. In section II, we describe our simulation protocol. Results gathered in section III are presented and discussed according to the following outline: in section III.1, the mechanical equilibrium achieved by the packings are carefully examined, as well as their homogeneity to ensure the absence of particle segregation; then section III.2 reports on the calculated solid fraction values and their variations as a function of particle bidispersity, angularity and aspect ratio; finally, section III.3 investigates the microstructure of our packings.

II Simulation protocol

The simulated systems are dense assemblies of $3000$ rigid frictionless particles of identical mass density $\rho$, interacting with each other through totally inelastic collisions. Particles are either spheres or $pinacoids$, the latter referring to a variety of convex polyhedra comprised of eight vertices, fourteen edges and eight faces as shown on figure 1. According to an extensive experimental study with various rock types mentioned in ref. Tourenq and Denis (1982), the pinacoid gives the best fit among simple geometries for an aggregate grain. As explained elsewhere Azéma et al. (2012); Camenen et al. (2012), a pinacoid has three planes of symmetry and is determined by four parameters, length $L$, width $G$, height $E$ ($L\geq G\geq E$) and angle $\alpha$ set to $60^{\circ}$ here, so that its volume $V$ is given by:

[TABLE]

Alternatively, a pinacoid may be determined by its size $d=\sqrt{G^{2}+L^{2}}$, corresponding to the diameter of its circumscribed sphere, supplemented with its aspect ratios $L/G$ and $G/E$. Figure 2 depicts pinacoids with various aspect ratios and Table 1 summarizes the size and aspect ratio of any particle used in the present study.

Binary mixes in which various proportions of small spheres have been substituted for large ones, or small isometric pinacoids for large isometric pinacoids, have been prepared to study the influence of particle size distribution and angularity on packing properties. Additionally, binary mixes comprising various proportions of flat, elongated, or flat & elongated pinacoids substituted for large isometric ones have been prepared to investigate the role played by particle aspect ratio. Table 2 summarizes the mean proportion by volume of small, elongated, flat or flat & elongated particles in each binary mix simulated.

For each binary mix, three cuboidal samples have been prepared following a pluviation protocol inspired by Laniel (2007) and described in details in Azéma et al. (2012); Camenen et al. (2012): spherical shells, each circumscribed to a randomly oriented particle, are first randomly dropped inside a vertical parallelepiped container and subsequently moved to the closest local minimum of potential energy; second, the spherical shells are removed, bi-periodic boundary conditions are substituted for the container vertical walls, and the same gravity $\vec{g}$ ([math],[math],$-g$) is applied to all particles until they find an equilibrium position under their own weight (examples of such an equilibrated packing is shown on Fig. 3). This protocol was selected for the ability of pluviation to achieve $RCP$ states of hard spheres with solid fraction repeatedly peaking at $0.64$ Bernal and Mason (1960); Macrae and Gray (1961); Scott and Kilgour (1969); Zhang et al. (2001); Silbert et al. (2002); Emam et al. (2005). Similarly, nonspherical frictionless particles subjected to the selected pluviation protocol were expected to achieve repeatedly states of maximum density characterized by a unique solid fraction value, provided that they are homogeneously spread and achieve a stable equilibrium without crystallization or segregation Roux (2004); Agnolin and Roux (2007).

All simulations were performed using the Non Smooth Contact Dynamic method (NSCD) Radjaï and Dubois (2011); Radjaï and Richefeu (2009); Jean (1999); Moreau (1994). This distinct element method (DEM), implemented in the LMGC90 sofware platform Dubois and Jean (2003), was successfully applied to a number of physical problems ranging from dense inertial flows Lois et al. (2005); Chevoir et al. (2001); Azéma et al. (2012) to quasistatic deformable packings Azéma et al. (2007); Azema et al. (2009); Estrada et al. (2008); Camenen et al. (2012). Basically, the equations of motion of a collection of rigid particles interacting through unilateral contacts with dry friction are first integrated over one time step. Hence, instead of accelerations and forces, the unknowns are particle velocities and percussions (also called $impulses$). The advantage of this integral formulation is that collision (shock) and lasting contact situations do not need to be distinguished any more, both giving rise to a percussion. Percussions, which may be defined as the integral of a force over one time step, are parameterized using normal ($e_{N}$) and tangential ($e_{T}$) restitution coefficients as well as a sliding friction coefficient ($\mu$). Second, at each time step, a geometrical contact detection algorithm identifies any potential contact situation—mainly vertex-face, edge-edge, edge-face or face-face contacts between convex polyhedra and a point between spheres—and describes it in terms of location and normal unit vector. Last, the equations of motion are solved at each time step by an iterative process using a non-linear Gauss-Seidel like method. This resolution is performed with respect to $complementarity$ $relations$ Moreau (1988) between relative velocities and percussions substituted for classical non-interpenetration (Signorini) and contact friction (Coulomb) constraints.

Simulation parameters are summarized in table 3.

III Results

We first check the quality of mechanical equilibrium achieved by our simulated packings, then we examine the variations of packing solid fraction as a function of particle size and shape. Finally, we investigate microstructural properties of our pinacoid packings to help understanding solid fraction variations.

III.1 Mechanical equilibrium and homogeneity

A packing of rigid particles achieves a mechanically stable equilibrium under the following conditions: $1$) the net force and net torque applied to each particle as well as its kinetic energy are negligible, and $2$) the potential energy of the collection of particles is minimal. For rigid spheres, all these conditions are met when the following criteria are fulfilled Camenen et al. (2012); Agnolin and Roux (2007):

[TABLE]

where $\sum F$, $\sum M$ and $E_{c}$ are respectively the net force, net torque and total kinetic energy of a particle, $P$ denotes the local average stress in its neighborhood and $F_{S}$ the net force exerted at contacts between particles and the bottom wall. Though our packings of frictionless spheres were found to meet these criteria, moderate deviations were observed for pinacoids. Indeed, the net force was found in the order of $10^{-3}d^{2}P$, whereas the net torque and kinetic energy were found in the order of $10^{-3}d^{3}P$ and $10^{-7}d^{3}P$ respectively. The fact that these results are slightly higher than previous ones Camenen et al. (2012) may be explained by shorter simulation durations (about $25s$ compared to about $35s$ in Camenen et al. (2012)), but this does not call into question the quality of mechanical equilibrium achieved (this was further checked upon continuing four simulations up to $35s$ - results are not shown here). Besides, all packings meet the requirement stated by equation (5) since the ratio of the total mass of particles to the net force $F_{S}$ exerted on the bottom wall varies between $0.9957$ and $1.0048$.

Next, we verify that no significant interpenetration occurs at interparticle contacts to ensure the relevance of solid fraction calculations. For sphere packings, the maximum calculated interpenetration is lower than $10^{-3}d$. Similarly, we verify that the interpenetration routinely generated at polyhedra contacts when applying the NSCD method Saussine (2004); Camenen et al. (2012) remains small enough. For this purpose, we apply the virtual works principle to generate a vertical expansion of each packing that will reset the vertical component of all interpenetrations to zero:

[TABLE]

where $k$ is the coefficient of vertical expansion of the packing (taken proportional to distance $z$ between the center of gravity of each particle and the bottom wall), $F_{C}^{z}$ and $\delta_{C}^{z}$ are the vertical components of contact forces and interpenetrations respectively, and both sums are calculated two particle size above the layered bottom to get rid of wall effect Camenen et al. (2012); Donev et al. (2004b). The mean $k$ value calculated over all pinacoid packings is $0.80\%$ and the sandard deviation is $0.47\%$. Upon assuming a similar coefficient value along $x$ and $y$ axes, interpenetrations at contacts are found to yield no more than $2$ to $3\%$ overestimation in solid fraction calculations, which is consistent with references Saussine (2004); Camenen et al. (2012).

Finally, we checked the homogeneity of the distribution of each population of particles once our packings have reached a stable mechanical equilibrium. Figure 4 depicts the proportions of large $P_{d}$, small $P_{d/3}$ and flat & elongated $P_{PO}$ particles as a function of distance $z$ from the bottom wall of their respective packings. The proportion profiles are almost constant, with only small deviations close to the bottom wall of bidisperse packings. Hence, we can conclude that each population of spheres or pinacoids is reasonably homogeneously spread in its respective packing.

III.2 Solid fraction

In this subsection, we investigate the effect of particle size distribution, angularity, and aspect ratio on the maximum solid fraction value reached by the packing once a stable mechanical equilibrium is achieved. For each packing, the solid fraction is computed upon evenly slicing the packing horizontally and performing analytical calculation of the volume of each sphere or each pinacoid present in each slice of known volume. Naturally, the mean solid fraction is calculated in the bulk, that is away from bottom wall and free surface.

III.2.1 Effects of particle size distribution and angularity

For various proportions of small particles $X_{S}$, Fig. 5 depicts maximal solid fraction values $\phi$ achieved by bidisperse packings of (a) spherical particles and (b) pinacoids. For spherical particles, excellent agreement can be observed between our results and those calculated from molecular dynamics Roux et al. (2007) or sphere mapping Farr and Groot (2010) simulations. The corresponding curves peak for $30\%$ of small spherical particles by volume, and their respective maximum values are $\phi=0.726\pm 0.004$ for the present study, $\phi=0.7207\pm 0.0004$ for ref. Roux et al. (2007) and $\phi=0.7324$ for ref. Farr and Groot (2010). A good agreement is also observed with experimental results reported in ref. Sedran et al. (1994) from bidisperse packings of ($d$, $d/4$) rounded aggregates densified inside a rigid cylinder through vertical taps and upper surface loading, though these are finite size packings of non-strictly bidisperse spherical particles. By contrast, similar experiments with bidisperse packings of ($d$, $d/2$) rounded aggregates yielded significantly lower maximum solid fractions, which may be explained by increased crowding of the local arrangement of large particles by small ones when their size ratio comes closer to $1$ Stovall et al. (1986).

For pinacoids, similar layout of the solid fraction curve is observed (Fig. 5b), with monodisperse isometric pinacoid packings ($X_{S}=0\%$ and $X_{S}=100\%$) achieving the minimum $\phi=0.676\pm 0.004$ and a peak $\phi=0.769\pm 0.001$ being reached when $X_{S}=30\%$. Interestingly, observe that this peak is significantly higher than the one calculated for sphere packings, suggesting that angular particles pack denser than spherical ones. This conclusion is consistent with previous experimentations reported by ref. Jaoshvili et al. (2010), showing that monodisperse plastic tetrahedra pack with a maximum solid fraction of $\phi=0.76\pm 0.02$ by extrapolation to infinite packing size, and with subsequent numerical work on Platonic solids reported in ref. Jiao and Torquato (2011). Note however that other previous works support experimentally Baker and Kudrolli (2010) or numerically Smith et al. (2010) the opposite conclusion illustrated on Fig. 5 by experimental results from ref. Sedran et al. (1994), according to which rounded particles would pack denser than angular ones. In fact, low experimental solid fraction values result from difficulties in overcoming interparticle friction while compacting a granular assembly Silbert et al. (2002); Agnolin and Roux (2007), which do not affect our frictionless pinacoids.

III.2.2 Shape effect

Figure 6 depicts maximum solid fraction $\phi$ versus the proportion by volume of $(a)$ elongated, $(b)$ flat and $(c)$ flat & elongated particles in mixtures with large isometric pinacoids. This figure shows that partial substitution of flat pinacoids for isometric ones results in packing solid fraction decrease by maximum $3\%$ for elongated and $8\%$ for flat & elongated particles. In both cases, the minimum solid fraction is achieved when the proportion of flat or flat & elongated particles reaches approximately $50\%$. Increasing further the proportion of flat particles raises the solid fraction value up to $\phi=0.676\pm 0.001$ for $100\%$ flat pinacoids and slightly above for $100\%$ flat & elongated pinacoids. By comparison, random jammed assemblies of $100\%$ oblate ellipsoids with an aspect ratio of $3$ were reported in ref. Donev et al. (2004b) to pack with a maximum solid fraction $\phi=0.67$. Contrarywise, partial substitution of elongated pinacoids for isometric ones does not seem to significantly affect the packing solid fraction which may be estimated at $\phi=0.683\pm 0.012$ for $100\%$ elongated pinacoids. Note that this estimate compares well with the $\phi=0.68$ solid fraction reported in ref. Kyrylyuk and Philipse (2011) for packings of $100\%$ spherocylinders with the same length to diameter aspect ratio of $2$. Upon applying the virtual works principle as depicted in equation (6), lower bounds accounting for particle maximum interpenetration may be calculated for the solid fraction of packings made of $100\%$ of the various particle shapes studied. Table 4 gathers these lower bounds and evidences that, regardless of their aspect ratios, pinacoids pack denser than spheres in the RCP state.

III.3 Microstructure

To shed some light on maximum solid fraction variations as a function of particle size, angularity and shape, we now investigate the microstructural properties of our packings. We first examine the existence of particle arrangement and then we explore the contact network.

III.3.1 Particle arrangement

The existence of translational arrangement is studied by means of the pair correlation function $g(r)$, which is calculated in the packing bulk (two particle size above the layered bottom to get rid of wall effect Camenen et al. (2012); Donev et al. (2004b)) using the expression detailed in ref. Allen and Tildesley (1989) (page $55$). Figure 7 depicts the variations of $g(r)$ in (a) bidisperse pinacoid packings and in (b) mixtures of isometric pinacoids with elongated, flat or flat & elongated particles. These curves clearly show that, for $r$ larger that $2R_{min}$ to $3R_{min}$, $g(r)$ no more significantly differs from $1$, meaning that the local distribution of particle centers is free of long-range translational ordering. Like tetrahedron packings Jaoshvili et al. (2010); Smith et al. (2010); Neudecker et al. (2013), pinacoid packings appear to be less correlated than sphere packings, because the former lack of rotational symmetry and angles between their adjacent flat faces induce frustration. Note that the peak visible on several curves for $r\simeq(2/\sqrt{6}+1/2)R_{min}\simeq 1.3R_{min}$ corresponds to the face-face contact between small particles as depicted on fig. 8.

Given their symmetry properties, pinacoids may potentially adopt the same orientation upon aligning one or more of their inertia axes, thus confering orientational order to the packing. To detect such an orientational order, the nematic order parameter $Q_{00}^{2}$ and the biaxial parameter $Q_{22}^{2}$ are computed in the packing bulk (see ref. Camenen et al. (2012); Camp and Allen (1997) for computation details). $Q_{00}^{2}$ assesses the highest level of alignement of a given inertia axis between all particles, either $\roarrow{u}$, $\roarrow{v}$, or $\roarrow{w}$ (see fig. 1), whereas $Q_{22}^{2}$ assesses overall alignment of particle inertia axes as a consequence of the pinacoid symmetry properties. Figure 9 depicts the variations of $Q_{00}^{2}$ for pinacoid packings incorporating various proportions by volume of small, elongated, flat or flat & elongated particles. It is remarkable that substituting small isometric pinacoids for large ones decreases the packing nematic order parameter whatever the substituted proportion. Similarly, substituting elongated particles for large isometric ones does not significantly increase the packing nematic order parameter. By contrast, substituting more than $30\%$ of flat or flat & elongated particles for large isometric ones increases the packing nematic order parameter beyond the isotropic-nematic transition Allen et al. (1993); Camp and Allen (1997), since flat particles tend to align their $\roarrow{w}$ inertia axis with the $\roarrow{z}$ direction as shown on fig. 10. This alignment evidences the reorientation of flat and flat & elongated pinacoids during the densification phase to form a layered structure which minimizes the potential energy of the equilibrated packing. This reorientation is facilitated by frustration release resulting from one particle dimension being significantly smaller than the other two. Besides, with the biaxial parameter never exceeding $0.05$ whatever the tested mixture, it should be pointed out that none of the tested mixture has reached the nematic-biaxial transition Camp and Allen (1997), in other words particle orientations remained random in the horizontal plane.

III.3.2 Contact network in the samples

In this subsection, we first check whether the contact network of frictionless isometric pinacoid packings fulfills the so-called isostatic conjecture. Next, we investigate whether gradually substituting small, elongated, flat, or flat & elongated pinacoids for isometric ones modifies the typology of the contact network.

The isostatic conjecture states that amorphous rigid packings of hard frictionless particles are isostatic Alexander (1998). In this statement, rigid packing refers to particle assembly which cannot be deformed subject to arbitrarily low stress without breaking any interparticle contact or deforming any particle, in other words packing in a mechanically stable equilibrium state. Furthermore, isostatic means that the total number of interparticle constraints equals the sum of particle degrees of freedom (DOF), that is the coordination number equals twice the particle DOF. Though disputed in the recent years Smith et al. (2010), authors have established the validity of this conjecture for random close packings of spheres Silbert et al. (2002); Agnolin and Roux (2007) and nontiling platonic solids such as tetrahedra, octaedra, icosaedra and dodecaedra Jaoshvili et al. (2010); Jiao and Torquato (2011), whereas packings of rounded particle shapes such as moderately oblate or prolate ellipsoids were found hypostatic Donev et al. (2004b).

The coordination number calculated here for monodisperse isometric pinacoid packings is $8.4\pm 0.2$. Note that this value lies between those observed for similar particle shapes, such as maximally randomly jammed monodisperse tetrahedra for which $N\approx 8.5$ and $N=8.6\pm 0.1$ are respectively reported in ref. Neudecker et al. (2013); Smith et al. (2010), or monodisperse octahedra for which $N\approx 7.7$ and $N=7.8\pm 0.1$ are respectively reported in ref. Smith et al. (2010); Jiao and Torquato (2011). However, this value is well below the isostatic number $N=12$ since it amalgamates diverse contact types corresponding to different numbers of constaints. Indeed, upon assigning $1$ constraint to each vertex-face or edge-edge contact (simple contact), $2$ constraints to each edge-face contact (double contact) and $3$ constraints to each face-face contact (triple contact) Jaoshvili et al. (2010); Jiao and Torquato (2011); Camenen et al. (2012), the number of constraints writes $N_{c}=N_{s}+2.N_{d}+3.N_{t}$, where $N_{s}$, $N_{d}$ and $N_{t}$ are, respectively, the numbers of simple, double and triple contacts. With $N_{s}=5.57\pm 0.06$, $N_{d}=2.12\pm 0.04$ and $N_{t}=0.72\pm 0.07$, we obtain $N_{c}=11.99\pm 0.35$. Despite the finite size of our packing periodic cell, observe that this value compares well with the isostatic number $N=12$ valid for infinite size packings. As a consequence, our monodisperse isometric pinacoid packings are reasonably isostatic as conjectured in reference Alexander (1998). Yet, it should be observed that the number of face-face contacts per particle ($N_{t}=0.72\pm 0.07$) is lower than those tabled in ref. Jiao and Torquato (2011) for tetrahedra ($2.21\pm 0.01$), even for octaedra ($1.44\pm 0.01$). As observed in ref Neudecker et al. (2013), we argue that the probability of perfect face-face alignment is low compared to that of either slightly shifted face-face or low angle edge-face contact. This is consistent with the finite slope/moderate first peak value of our pair correlation curves (see Fig. 7) and our significantly higher number of edge-face contacts per particle ($N_{d}=2.12\pm 0.04$) compared to those reported for tetrahedra ($0.98\pm 0.01$) and octaedra ($1.38\pm 0.01$) in reference Jiao and Torquato (2011).

To get insight into how the contact network fluctuates with increasing proportions of small or non-isometric particles, it is worth examining Fig. 11 depicting the variations of the coordination numbers of all, simple, double and triple contacts as a function of the proportion by volume of small, elongated, flat, or flat & elongated pinacoids. When gradually substituting small isometric pinacoids for large ones, the coordination number first reaches a minimum $N=7.4\pm 0.2$ for $X_{S}=13\%$ and then gradually increases back to its monodisperse value, as shown on Fig. 11a. Interestingly, similar coordination number fluctuations may be observed with simple contacts (Fig. 11b), whereas the coordination numbers of double and triple contacts vary conversely (Fig. 11c and 11d). In fact, small pinacoids are trapped inside the excluded volume of large ones. Until their proportion is sufficient for completely filling this volume, which corresponds to the achievement of the packing maximum solid fraction ($X_{S}=30\%$, see Fig. 5b), their low steric hindrance allows them to rotate and establish more stable face-face or edge-face contacts with large ones. Contrarywise, gradually substituting elongated, flat, or flat & elongated pinacoids for isometric ones tends to increase the coordination number $N$ of our pinacoid packings up to, respectively, $8.9\pm 0.1$, $9.1\pm 0.2$, and $9.2\pm 0.2$. Actually, these non-isometric particles are too large to fit in the excluded volume of large isometric pinacoids while leaving there arrangement undisturbed. Besides, contrary to elongated or isometric pinacoids, the adjacent trapezoidal flat faces of flat and flat & elongated pinacoids do not intersect at right angle and these particles have a smaller thickness than the former. Stated otherwise, flat and flat & elongated pinacoids have a broken angular symmetry and reduced external surface (at least reduced triangular flat faces for flat & elongated pinacoids), causing the coordination number of face-face contacts to decrease to the benefit of simple contacts with increasing proportions of flat or flat & elongated pinacoids in the packing (see Fig. 11b and d). Furthermore, for increasing proportions of flat or flat & elongated pinacoids, observe that the coordination number of simple contacts increases quicker than the coordination number of face-face contacts decreases, which is a consequence of packing isostaticity : whatever the proportion of small, elongated, flat, or flat & elongated particles, all our packings remain isostatic with $12.02\pm 0.35$ constraints per particle.

Finally, Fig. 12 depicts the variations of the coordination number as a function of the nematic order parameter. This figure suggests that the coordination number increases continuously with orientational ordering in the packings and irrespective of particle size or aspect ratio. An exponential fit of the following form was adjusted to the point cloud:

[TABLE]

where [$N(0),N(1),Q^{2}_{00,c}$] were calculated respectively as ($7.4,9.2,0.2$) with $R^{2}=0.92$. This fit suggests that orientationnally ordered pinacoid packings have coordination numbers in excess of about $9$ and, conversely, that packings may be considered orientationnally disordered for coordination numbers below roughly $8.5$, corresponding to $Q^{2}_{00}\leq 0.2$.

IV Conclusion

The properties of dense packings of spheres or pinacoids compacted under their own weight have been investigated using three-dimensional Non Smooth Contact Dynamic simulations. Various proportions by volume of small, elongated, flat, or flat & elongated pinacoids were substituted for large isometric ones in order to understand how polydispersity and shape affect their solid fraction and microstructural properties. Numerical simulations show that disordered assemblies of frictionless pinacoids, were they monodisperse or bidisperse, pack with a higher solid fraction than corresponding assemblies of spherical or rounded particles, thus fulfilling the analogue of Ulam’s conjecture for random packings proposed in ref. Jiao and Torquato (2011). This seeming discrepancy with experimental results reported in ref. de Larrard (1999); Baker and Kudrolli (2010) is believed to lie with difficulties in overcoming interparticle friction through experimental densification processes. Moreover, solid fraction increases further with bidisperse particles and peaks when the proportion of small ones reaches $30\%$, achieving $\phi=0.726\pm 0.004$ and $\phi=0.769\pm 0.001$ respectively for spheres and pinacoids. Contrarywise, partial substitution of flat pinacoids for isometric ones results in packing solid fraction decrease by a maximum $8\%$, especially when flat particles are also elongated. Minimum solid fraction is achieved when the proportion by volume of flat or flat & elongated pinacoids reaches $50\%$. Nevertheless, particle shape seems to play a minor role on packing solid fraction compared to polydispersity. Additional investigations focused on packing microstructure confirm that, with $12.02\pm 0.35$ constraints per particle, pinacoid packings fulfill the isostatic conjecture and that they are free of order except beyond $30$ to $50\%$ by volume of flat or flat & elongated polyhedra in the packing. This order increase progressively takes the form of a nematic phase as flat or flat & elongated particles reorientate so that their largest projected area is horizontal to minimize the packing potential energy. Simultaneously, this reorientation seems to increase the solid fraction value slightly above the maximum achieved by monodisperse isometric pinacoids, as well as the coordination number. Finally, partial substitution of elongated pinacoids for isometric ones has limited effect on packing solid fraction or order.

Acknowledgements.

The authors thank the team running the Centre de Calculs Intensifs des Pays de la Loire (CCIPL) for providing the calculation ressources through the MTEEGD project. They are also grateful to their colleagues J.N. Roux, N. Roquet and P. Richard for helpful conversations. Many thanks to the LMGC90 team in Montpellier for making their software plateform freely available.