Random close packing of disks and spheres in confined geometries

PHYSICAL REVIEW E 80, 051305 共2009兲

Random close packing of disks and spheres in confined geometries Kenneth W. Desmond and Eric R. Weeks Department of Physics, Emory University, Atlanta, Georgia 30322, USA 共Received 3 March 2009; revised manuscript received 22 July 2009; published 30 November 2009兲 Studies of random close packing of spheres have advanced our knowledge about the structure of systems such as liquids, glasses, emulsions, granular media, and amorphous solids. In confined geometries, the structural properties of random-packed systems will change. To understand these changes, we study random close packing in finite-sized confined systems, in both two and three dimensions. Each packing consists of a 50-50 binary mixture with particle size ratio of 1.4. The presence of confining walls significantly lowers the overall maximum area fraction 共or volume fraction in three dimensions兲. A simple model is presented, which quantifies the reduction in packing due to wall-induced structure. This wall-induced structure decays rapidly away from the wall, with characteristic length scales comparable to the small particle diameter. DOI: 10.1103/PhysRevE.80.051305

PACS number共s兲: 45.70.⫺n, 61.20.⫺p, 64.70.kj

I. INTRODUCTION

Random close packing 共rcp兲 has received considerable scientific interest for nearly a century dating back to the work of Westman in 1930 关1–6兴 primarily due to the relevance rcp has to a wide range of problems, including the structure of living cells 关7兴, liquids 关8,9兴, granular media 关10–13兴, emulsions 关14兴, glasses 关15兴, amorphous solids 关16兴, jamming 关17兴, and the processing of ceramic materials 关18兴. Typically, one defines rcp as a collection of particles randomly packed into the densest possible configuration. More rigorous definitions are available 关7兴, but it is generally accepted that the rcp density of a packing of spheres is ␾rcp ⬇ 0.64. Packings can have other rcp densities when the particles are polydisperse mixture of spheres 关19–24兴, nonspherical in shape 关25–29兴, or confined within a container that is comparable in size to a characteristic particle size 关18,30–39兴. While most studies of rcp focus on infinite systems, real systems have boundaries and often these boundaries are important as highlighted by Carman in 1937 关30兴. In the experiments by Carman, the packing fraction dependence on container size was measured for spheres poured into a cylindrical container and shaken for sufficiently long enough time to reach a very dense state. It was found that the packing fraction decreases with container size, which was attributed to the boundaries, altering the structure of the packing in the vicinity of the wall. Since the work of Carman, there have been many other studies, which have investigated rcp in confined systems 关6,35,40,41兴. These studies have shown that near the boundary, particles tend to pack into layers, giving rise to a fluctuating local porosity with distance from the wall, ultimately affecting the macroscopic properties of highly confined systems. Other studies have examined the packing of granular particles in narrow silos, focusing on the influence of confinement on stresses between particles and the wall 关42–45兴. Nearly all of these studies did not directly measure the local packing or any local packing parameters with relation to distance from the side wall, with the exception of a few experiments that used x-ray imaging to view the structure of confined packings. In these experiments, the packings were 1539-3755/2009/80共5兲/051305共11兲

monodisperse, facilitating highly ordered packing near the boundary, with measurements carried out at only a few different container size to particle size ratios 关37,41兴. Even with the history of work on the study of rcp in confined geometries, there is little known about how sensitive the structure of the packing near the boundary is to small changes in the confining width. For example, prior work found nonmonotonic dependence of ␾rcp on container size but only at extremely small containers with narrow dimensions h only slightly larger than the particle diameter d, that is, h ⬇ 3d or smaller 关30,36,38兴. However, their data were not strong enough to look for such effects at larger container sizes. Additionally, primarily only confined monodisperse systems have received much attention, and these systems are susceptible to crystallization near flat walls, which greatly modify the behavior 关46兴. 共One group did study binary systems, but they were unable to directly observe the structure 关35兴.兲 Furthermore, two-dimensional 共2D兲 confined systems have not been studied systematically, although they are relevant for a wide range of granular experiments 关47兴. In this paper, we address these questions using computersimulated rcp packings in confined geometries. In particular, we study binary mixtures to prevent wall-induced crystallization 关48–50兴. We create 2D and three-dimensional 共3D兲 packings with flat confining walls. In some cases, the system is confined only along one dimension 共with periodic boundaries in the other directions兲, and in other cases we confine the sample along all directions. Our simulations are carried out at many different and very closely spaced confining thicknesses, spanning a large range of values to elucidate the effects small changes in confining thickness has on the structure. We find that confinement significantly modifies the rcp states, with lowered values for ␾rcp, reflecting an inefficient packing near the walls. This inefficient packing persists several particle diameters away from the wall, although its dominant effects are only within 1 to 2 diameters. The behavior of ␾rcp is not monotonic with increasing sample thickness, reflecting the presence of boundary layers near the walls. Understanding the character of random close packing in confined geometries may be relevant for non-close-packed confined situations 关51兴. For example, when a liquid is con-

051305-1

©2009 The American Physical Society

PHYSICAL REVIEW E 80, 051305 共2009兲

KENNETH W. DESMOND AND ERIC R. WEEKS

fined, its structure is dramatically changed; particles form layers near the wall, which ultimately affects the properties of the liquid 关52–57兴. The shearing of confined dense colloidal suspensions shows the emergence of new structures not seen before 关58兴. The flow of granular media through hoppers 关59,60兴 or suspensions through constricted microfluidic and nanofluidic devices 关61–64兴 can jam and clog, costing time and money. One of our own motivations for this work was to help us understand prior experiments by our group, which studied the confinement of colloidal particles 关65兴. A dense suspension of colloidal particles behaves similarly to a glass 关66兴. For traditional glass formers, many experiments have studied how confinement modifies the glass transition; samples, which have a well-characterized glass transition in large samples, show markedly different properties when confined to small samples 关50,51,67–72兴. In our experimental work, the colloidal particles had much slower diffusion rates when confined between two parallel walls 关65兴. However, the experiments were difficult and we only examined behavior of a dense suspension at a few specific thicknesses. As noted above, in this current work we investigate how particles pack for a finely spaced set of thicknesses, to look for nonmonotonic behavior of the packing that might have been missed in the experiment. A second related question is whether confinement effects on glassy behavior are due to boundary effects or finite-size effects 关73兴. Our results show that boundaries significantly modify the packing, which may in turn modify the behavior of these confined molecular systems 关51兴. The paper is organized as follows. Section II outlines the algorithm we use to generate confined rcp states. Section III shows how the total packing fraction, particle number density, and local order of confined rcp states change with confining thickness and distance from the confining boundary. Finally, Sec. IV provides a simple model that predicts the packing fraction dependence with confinement. II. METHOD

Our aim is to quantify how a confining boundary alters the structure of rcp disks in 2D and spheres in 3D and, in particular, to study how this depends on the narrowest dimension. This section presents our algorithm for 2D packings first and then briefly discusses differences for the 3D algorithm. In 2D, our system consists of a binary mixture of disks containing an equal number N / 2 large disks of diameter dl and small disks of diameter ds with size ratio ␴ = dl / ds = 1.4. For each configuration, disks are packed into a box of dimensions Lx by Ly. For most simulations we discuss, there is a periodic boundary condition along the x direction and two fixed hard boundaries 共walls兲 along the y direction; although as discussed below, in some cases we consider periodic boundaries in all directions or fixed boundaries in all directions. Each configuration is generated using a method adapted from Xu et al. 关74兴, which is an extension of a method proposed by Clarke and Wiley 关75兴. This method is briefly sum-

FIG. 1. A flow chart outlining our algorithm for computing rcp configurations.

marized in Fig. 1. Infinitesimal particles are placed at random 关76兴 in the system, gradually expanded and moved at each step to prevent particles from overlapping. When a final state is found such that particles can no longer be expanded without necessitating overlap, the algorithm terminates. Near the conclusion of the algorithm, we alternate between expansion and contraction steps to accurately determine the state. In particular, while the final state found is consistent with hard particles 共no overlaps allowed兲, the algorithm uses a soft potential at intermediate steps 关74兴, given by

⑀ V共rij兲 = 共1 − rij/dij兲2⌰共1 − rij/dij兲, 2

共1兲

where rij is the center to center distance between two disk i and j, ⑀ is a characteristic energy scale 共⑀ = 1 for our simulations兲, dij = 共di + d j兲 / 2, and ⌰共1 − rij / dij兲 is the Heaviside function, making V nonzero for rij ⬍ dij. Simulations begin by randomly placing disks within a box of desired dimensions and boundary conditions with the initial diameters chosen such that ␾initial Ⰶ ␾rcp. In the initial state, particles do not overlap and the total energy E = 0. Next, all disk diameters are slowly expanded subject to the fixed size ratio ␴ = 1.4 and ␾ changing by ␦␾ per iteration; we start with ␦␾ = 10−3. After each expansion step, we check if any disks overlap by checking the condition 1 − rij / dij ⬎ ⑀r = 10−5 for each particle pair. Below this limit, we assume the overlap is negligible. If any particles do overlap 共E ⬎ 0兲, we use the nonlinear conjugate gradient method 关77兴 to decrease the total energy by adjusting the position of disks so they no longer overlap 共E = 0兲. In practice, one energy minimization step does not guarantee that we have reached a minimum within the desired numerical precision. Thus, this step can be repeated to further reduce the energy if E ⬎ 0. We judge that we have reached a nonzero local mini-

051305-2

PHYSICAL REVIEW E 80, 051305 共2009兲

RANDOM CLOSE PACKING OF DISKS AND SPHERES IN…

FIG. 2. 共Color online兲 Illustrations of 2D and 3D configurations generated using the algorithm described in Sec. II. 共a兲 2D configuration for h = 10. 共b兲 2D configuration with h = 20. 共c兲 3D configuration with h = 5, where blue 共dark gray兲 represents big particles and green 共light gray兲 represents small particles. 共d兲 3D configuration with h = 10.

mum if the condition 储ⵜE储 / 共2N兲 ⬍ ⑀E = 10−7 is found, where 储ⵜE储 is the magnitude of the gradient of E. Physically speaking, this is the average force per particle, and the threshold value 共10−7兲 leads to consistent results. If we have such a state with E ⬎ 0, this is not an rcp state as particles overlap. Thus, we switch and now slowly contract the particles until we find a state where particles again no longer overlap 共within the allowed tolerance兲. At that point, we once again begin expansion. Each time we switch between expansion and contraction, we decrease ␦␾ by a factor of 2. Thus, these alternating cycles allow us to find an rcp state of nonoverlapping particles 共within the specified tolerance兲 and determine ␾rcp to high accuracy. We terminate our algorithm when ␦␾ ⬍ ␦␾min = 10−6. In practice, we have tested a variety of values for the thresholds ⑀r, ⑀E, and ␦␾min and find that our values guarantee reproducible results as well as reasonably fast computations. Our algorithm gives an average packing fraction of ␾rcp = 0.8420⫾ 0.0005 for 40 simulated rcp states, containing 10 000 particles with periodic boundary conditions along both directions. Our value of ␾rcp is in agreement with that found by Xu et al. 关74兴. The above procedure is essentially the same as Ref. 关74兴; we modify this to include the influence of the boundaries. To add in the wall, we create image particles reflected about the position of the wall; thus, particles interact with the wall using the same potential 关Eq. 共1兲兴. Additionally, we wish to generate packings with prespecified values for the final confining height h = Ly / ds. 共This allows us to create multiple rcp configurations with the same h.兲 We impose h by affinely scaling the system after each step, so that the upper boundary is adjusted by Ly = hds and each disk’s y coordinate is multiplied by the ratio Ly,i+1 / Ly,i, where Ly,i and Ly,i+1 are the confining widths between two consecutive iterations. Thus, while ds gradually increases over the course of the simulation, Ly increases proportionally so that the nondimensional ratio h is specified and constant. Some examples of our final rcp states are shown in Fig. 2. To ensure we will have no finite-size effects in the periodic direction, we examined ␾rcp for different h and Lx and found ␾rcp共h兲 to be independent of Lx for 3 ⱕ h ⱕ 30 if Lx / ds ⬎ 40. Thus, we have chosen N for each simulation to guarantee Lx / ds ⬇ 50. In 3D, our system consists of a binary mixture of spheres, containing an equal number N / 2 large spheres of diameter

dl and small spheres of diameter ds with a size ratio ␴ = dl / ds = 1.4. Spheres are packed into a box of dimensions Lx by Ly by Lz, with periodic boundaries along the x and z directions and a fixed hard boundary along the y direction. Each configuration is generated using the same particle expansion and contraction method described above and the same initial values for ␦␾ and the terminating conditions. For each configuration Lx = Lz, h = Ly / ds, and N is chosen so that Lx / ds ⬎ 10. Our choice of Lx / ds ⬎ 10 is not large enough to avoid finite effects. However, in order to acquire the large amount of data needed in a reasonable amount of time, we intentionally choose a value of Lx / ds below the finite-size threshold. Trends observed in the 2D analysis will be used to support that any similar trends seen in 3D are real and not the result of the finite periodic dimensions. Note that in 3D we will show cases where h ⬎ Lx / ds, resulting in the confining direction being larger than the periodic direction, and this may affect the structure of final configurations; however, we will not draw significant conclusions from those data. Overall, it is not known if this algorithm produces mathematically rigorously defined random close-packed states 关7,17,74,78兴. However, the goal of this paper is to determine empirically the properties of close-packed states in confinement, and we are not attempting to extract mathematically rigorous results. For example, we are not as interested in the specific numerical values of ␾rcp that we obtain, but rather the qualitative dependence on h. As noted in the introduction, different computational and experimental methods for creating rcp systems have different outcomes, and so it is our qualitative results we expect will have the most relevance. Note that for the remainder of this paper, we will drop the subscript rcp, and it should be understood that discussions of ␾ refer to the final state found in each simulation run ␾rcp共h兲. III. RESULTS A. 2D systems

We begin by generating many 2D configurations with h between 3-30 and computing the packing fraction for each, as shown by the black curve in Fig. 3. This plot shows that confinement lowers ␾, with the influence of the walls being increasingly important at lower h. The lowering of ␾ with confinement is most likely due to structural changes in the packing near the confining boundary. We know that any alteration in particle structure from a rcp state must be “near” the wall because as h → ⬁, we expect to recover a packing fraction of ␾rcp, implying that in the infinite system the “middle” of the sample is composed of an rcp region. Extrapolating the data in Fig. 3 to h → ⬁, we find ␾h→⬁ = ␾rcp = 0.842, which is essentially a test of our method. The extrapolation 共red curve in Fig. 3兲 was carried out by assuming that to first order ␾ ⬃ ␾h→⬁ − C / h for large h, where ␾h→⬁ = ␾rcp 共the bulk value for the rcp packing兲 and C is a fitting parameter. The data in Fig. 3 begin to deviate from the fit for h ⱗ 6 and, furthermore, ␾共h兲 is not monotonic. While some of the variability is simply noise due to the finite number of disks N used in each simulation, some of the variability is real. The inset in Fig. 3 shows a magnified view of the region

051305-3

PHYSICAL REVIEW E 80, 051305 共2009兲

KENNETH W. DESMOND AND ERIC R. WEEKS

FIG. 3. 共Color online兲 The black curve is the average packing fraction ␾ found by averaging at least ten 2D configurations together for various confining widths h; recall that h has been nondimensionalized by ds, the small particle diameter. The smooth red curve 共dark gray兲 is a fit using Eq. 共5兲, which finds ␾rcp = 0.842 in the limit h → ⬁; the value for ␾rcp is indicated by the black dashed line. The green 共light gray兲 data points are ␾共h兲 computed for many configurations with the confining wall replaced by a periodic boundary. The inset is a magnified view of the region for h ⱕ 6 to better show the large variations within this range. The vertical lines in the inset are located at “special” h values, where peaks and plateaus appear.

3 ⱕ h ⱕ 6. The vertical lines in this inset are located at specific values of h that can be expressed as the integer sums of the two-particle diameters. For instance, the first vertical line near the y axis is located at h = 2ds + dl. These lines are placed at some h values, where ␾共h兲 has notable spikes or plateaus. These lines suggest that there exist special values of h, where the confining thickness is the right width so that particles can pack either much more efficiently or much less efficiently than nearby values of h. Intriguingly, these special h values do not appear to be as pronounced at all possible integer sums, but instead only the selected few are drawn in the figure. However, given the apparent noisy fluctuations 共despite averaging over a very large number of simulations兲, we cannot completely rule out that local maxima and minima might also exist at other combinations of ds and dl. Somewhat surprisingly, we do not observe large peaks corresponding to integer combinations of 共冑3 / 2兲ds or 共冑3 / s兲dl, which would suggest hexagonal packing, the easiest packing of monodisperse disks in 2D; whereas, the observed peaks of ␾共h兲 suggest squarelike packing. To measure structural changes in particle packing as a result of confinement, we start by examining the variations in the local number density ␳ with distance y from the confining wall. We define ␳ to be the average number of particles per unit of area along the unconfined direction. For a given location y, we count the number of particles in a region of area Lx⌬y and divide by this area, choosing ⌬y to be of a size such that the results do not depend sensitively on the choice but also so that we can get reasonably localized information. Figure 4 is a plot of ␳共y兲 for 100 configurations averaged together at h = 30. This plot shows oscillations in particle density, which decay to a plateau. The oscillations near the wall are indicative of particles layering in bands. Above y ⲏ 6ds, noise masks these oscillations. This supports our interpretation that confinement modifies the structure near the walls but not in the interior. Furthermore, the rapid-

FIG. 4. 共Color online兲 A plot of the number density ␳共y兲 for 100 2D configurations at h = 30 averaged together. The plot is constructed by treating the small and big particles separately and using bins along the confining direction of width ␦y = 0.1ds.

ity of the decay to the plateau seen in Fig. 4 suggests that confinement is only a slight perturbation to systems with overall size h ⲏ 6. The details of the density profiles in Fig. 4 also suggest how particles pack near the wall. The small particle density 共solid line兲 has an initial peak at y = 0.5ds, indicating many small particles in contact with the wall, as their centers are one radius away from y = 0. Likewise, the large particle density 共dashed line兲 has its initial peak at y = 0.7ds = 0.5dl, indicating that those particles are also in contact with the wall. This is consistent with the pictures shown in Figs. 2共a兲 and 2共b兲, where it is clear that particles pack closely against the walls. Examining again the small particle number density in Fig. 4 共solid line兲, the secondary peaks occur at y = 1.5ds and y = 1.9ds = 0.5ds + 1.0dl, which is to say either one small particle diameter or one large particle diameter further away from the first density peak at y = 0.5ds. This again is consistent with particles packing diameter to diameter rather than “nesting” into hexagonally packed regions. Similar results are seen for the large particles 共dashed line兲, which have secondary peaks at y = 1.0ds + 0.5dl and y = 2.1ds = 1.5dl. To confirm that these density profile results apply for a variety of thicknesses h, and more importantly to see how these results are modified for very small h, we use an image representation shown in Fig. 5. To create this image, density distributions of different h are each separately rescaled to a maximum value of 1. Every data point within each distribution is then made into a grayscale pixel indicating its relative value; black is a relative value of 1, and white is a relative value of 0. The vertical axis is the confining width and the horizontal axis is the distance y from the bottom wall. Each horizontal slice 共constant h兲 is essentially the same type of distribution shown in Fig. 4. The white space on the right side of Fig. 5 arises because the distribution is only plotted for the range 0 ⱕ y ⱕ h / 2. The distributions are symmetric about y = h / 2 and by averaging the distribution found for the range 0 ⱕ y ⱕ h / 2 with the distribution found for the range h / 2 ⱕ y ⱕ h, the statistics are doubled. The areas shown in the insets are magnified views, where the full range 0 ⱕ y ⱕ h is being shown. In Fig. 5, there are vertical strips of dark areas, once again indicating that particles are forming layers. The width of

051305-4

PHYSICAL REVIEW E 80, 051305 共2009兲

RANDOM CLOSE PACKING OF DISKS AND SPHERES IN…

␭=

FIG. 5. An image representation constructed for the purpose of comparing 2D ␳共y兲 distributions at many different h. The intensities have been logarithmically scaled. The vertical pixel width is 0.1 and for the left plot the horizontal pixel width is 0.2 and for the right plot the horizontal pixel width is 0.14.

these strips widens and the intensity lessens farther from the wall. In each plot, the first vertical black strip is sharply defined and located at one-particle radius, illustrating that small and big particles are in contact with the wall. Finally, the location and width of each layer remain essentially the same for different h, suggesting that layering arises from a constraint imposed by the closest boundary. Given that the first layer of particles always packs against the wall, this imposes a further constraint on how particles pack in the nearby vicinity. The consistency in the location and width of the second layer for all h demonstrates that the constraint of the first layer always produces a similar packing in the second layer, essentially independent of h. Continuing this argument, each layer imposes a weaker constraint on the formation of a successive layer, allowing for the local packing to approach rcp far from the wall. In the magnified views of Fig. 5, the vertical dark lines show the layering of particles induced by the left boundary and the angled dark lines show the layering of particles induced by the right boundary. We see that for small h, these sets of lines overlap and intersect, meaning that there is a strong influence from one boundary on the packing within the layers produced by the other boundary. This may explain the variations seen in ␾共h兲 for small h in Fig. 3. In particular, it is clear that at certain values of h, the layers due to one wall are coincident with the layers due to the other wall, and this suggests why ␾共h兲 has a higher value for that particular h. Given that the layer spacings correspond to integer combinations of ds and dl, the coincidence of layers from both walls will correspond to integer combinations of ds and dl, and this thus gives insight into the peak positions shown in the inset of Fig. 3. As described above, the influence of the walls diminishes rapidly with distance y away from the wall. In particular, for the local number density ␳共y兲, we observe that the asymptotic limit ␳共y → ⬁兲 = 0.362 for the curves shown in Fig. 4 is in agreement with the theoretical number density of an rcp configuration ␳rcp = 4␾rcp / ␲共1 + ␴兲. To quantify the approach to the asymptotic limit, we define a length scale from a spatially varying function f共y兲 using

冕 冕

y关f共y兲 − f共y → ⬁兲兴2dy .

共2兲

关f共y兲 − f共y → ⬁兲兴 dy 2

In this equation, f共y兲 is an arbitrary function, where the value of ␭ quantifies the weighting of f共y兲. For simple exponential decay f共y兲 = Ae共−y/␭⬘兲, Eq. 共2兲 gives ␭ = ␭⬘ / 2. Using f共y兲 = ␳共y兲, we find ␭ = 0.85ds and ␭ = 0.72ds for the small particle curve and big particle curve in Fig. 4, respectively, suggesting that the transition from wall-influenced behavior to bulk rcp packing happens extremely rapidly. To further investigate the convergence of the local packing to rcp more closely, we analyze the local bond order parameters ␺n, which for a disk with center of mass ri are defined as

␺n共ri兲 =

1 nb

冏兺 j

冏

eni␪共rij兲 .

共3兲

The sum is taken over all j particles that are neighbors of the ith particle, ␪共rij兲 is the angle between the bond connecting particles i and j and an arbitrary fixed reference axis, and nb is the total number of i-j bonds 关79兴. 共These are not physical bonds but indicate that two particles are nearest neighbors, where the definition of nearest neighbor is set by the first minimum of the pair-correlation function.兲 The magnitude of ␺2n is bounded between zero and one; the closer the magnitude of ␺2n is to 1, the closer the local arrangement of neighboring particles is to an ideal n-sided polygon. Figures 6共a兲–6共c兲 are drawings illustrating the concept of ␺2n using a 2D configuration with h = 10. Particles with larger ␺2n are drawn darker. These figures have no large clusters of dark colored particles, demonstrating that there are no large crystalline domains 共i.e., particles are randomly packed兲. For a highly ordered monodisperse packing, 具␺26典 would be the most appropriate choice for measuring order because of the ability for monodisperse packings to form hexagonal packing. However, for a bidisperse packing with size ratio ␴ = 1.4, the average number of neighbors a small particle will have is 5.5 and the average number of neighbors big particles will have is 6.5. Therefore, a bidisperse packing of this kind will have a propensity to form local pentagonal, hexagonal, and heptagonal packing, and to properly investigate how the local packing varies we examine 具␺25典, 具␺26典, and 具␺27典. We compute the average values 具␺25典, 具␺26典, and 具␺27典 for all configurations, as a function of y, and averaging together all 具␺2n典 distributions for configurations with h ⱖ 16 to improve statistics. This averaging can be justified by considering that oscillations in ␳共y兲 in Fig. 4 for y / ds ⬎ 10 are quite small. Thus, this averaging improves our statistics for the range 0 ⬍ y / ds ⬍ 5, where the largest oscillations occur, without skewing the data. In the end, nearly 10 000 configurations are averaged together, producing the curves shown in Figs. 7共b兲–7共d兲. This figure shows the spatial variations of 具␺25典, 具␺26典, and 具␺27典 for small and big particles separately and both particles combined. All curves show fluctuations that decay with distance from the wall and show local order within and

051305-5

PHYSICAL REVIEW E 80, 051305 共2009兲

KENNETH W. DESMOND AND ERIC R. WEEKS

2

2

2

2

2

2

FIG. 6. 共Color online兲 Drawings illustrating the conceptual meaning of 共a兲 ␺25, 共b兲 ␺26, and 共c兲 ␺27. Darker colored particles have neighbors that are packed more like an ideal regular n-sided polygon as compared to lighter drawn particles. The configuration of particles is the same for all panels and are drawn from a simulation with h = 10. Note in 共b兲 that there are no large patches of high ␺26, demonstrating that there are no large crystalline domains.

between layers. Figure 7共a兲 has been added so comparisons between the locations of the oscillations in ␳共y兲 and 具␺2n典共y兲 can be made. Each successive layer has less orientational order than the previous layer with 具␺2n典 eventually decaying to an asymptotic limit. To characterize a length scale for these curves, we compute ␭ using Eq. 共2兲 for each curve shown in Figs. 7共b兲–7共d兲. From the nine curves, we find that the mean value of ␭ = 共1.00⫾ 0.24兲ds. The length scales found for these curves are once again less than the largest particle diameter. No striking difference is found between the different order parameters or between the different particle sizes; specific values of ␭ are given in the figure caption. 共Note that the asymptotic limits of all 具␺2n典 plots are in agreement with the average values found for 40 unconfined 10 000 particle simulations averaged together, confirming that the local packing converges to an rcp arrangement far from the walls.兲 Next, we wish to distinguish the structural influence of the flat wall from the finite-size effects. We perform simulations, where the confining wall is replaced by a periodic boundary with periodicity h; thus, particles cannot form layers. In this case, the packing fraction still decreases as h is decreased, as shown by the green curve 共light gray兲 in Fig. 3, although the effect is less striking than for the case with walls 共black curve兲. A likely explanation for the decrease in ␾ with confinement is the long-range structural correlations imposed along the constricted direction; in other words, if there is a particle located at 共x , y兲 that particle is mirrored at 共x , y − h兲 and 共x , y + h兲 by the periodic boundary. We know from the pair-correlation function 关19,21兴 of rcp configurations that structural correlations exist over distances of many particle

FIG. 7. 共Color online兲 共a兲 A plot of the local number density ␳共y兲 for 2D configurations of big and small particles separately. 关共b兲–共d兲兴 Plots of 具␺2n典共y兲 for small 共green/light gray兲 and big particles 共blue/dark gray兲 separately and both sizes together 共light purple/medium gray兲 where 共b兲 is 具␺25典, 共c兲 is 具␺26典, and 共d兲 is 具␺27典. The length scales determined from these curves for small, large, and both species are ␭5,s = 1.2, ␭5,l = 1.1, ␭5,b = 1.4, ␭6,s = 0.8, ␭6,l = 0.9, ␭6,b = 0.8, ␭7,s = 1.1, ␭7,l = 0.7, and ␭7,b = 1.0 共all in terms of ds兲.

diameters, although of course these are weak at larger distances. Thus, the periodicity forces a deviation from the ideal rcp packing that becomes more significant as h decreases. By definition, rcp is the most random densely packed state and, thus, any perturbations away from this state must have a lower packing fractions. However, this is not nearly as significant as the constraint imposed by the flat wall, as is clear comparing the green 共light gray兲 data and the black data in Fig. 3. B. 3D systems

We next consider 3D confined systems. We start by investigating ␾共h兲 shown as the black points in Fig. 8. As observed in the 2D case, ␾ is reduced as a result of confine-

FIG. 8. 共Color online兲 The black data points are the average packing fractions of 3D configurations at various h. The red 共dark gray兲 curve is a fit to the model 共5兲. For each h, at least ten configurations were averaged together.

051305-6

PHYSICAL REVIEW E 80, 051305 共2009兲

RANDOM CLOSE PACKING OF DISKS AND SPHERES IN…

In 2D, we also noted that the structure is modified near the wall, as measured by the ␺2n order parameters. To investigate structural ordering in 3D, we use a local structural parameter sensitive to ordering 关80,81兴. We start by defining qˆi,6 =

FIG. 9. 共Color online兲 共a兲 A plot of ␳共y兲 for 3D configurations for small and big particles separately. The plot is constructed using bins of width ␦y = 0.1ds along the confining directions. 共b兲 A plot of the average number of ordered bonds 具Nb典共y兲.

ment. However, unlike the 2D system, there does not appear to be a series of “special values” of h that give rise to peaks and plateaus other than a hump near h = 3.75. The lack of substructure may be due to the smaller size in x and z, in contrast with the 2D simulations, which had large sizes in the unconfined direction. Next we investigate the local number density ␳共y兲 共the average number of particles per unit area along the unconfined directions兲 for h = 25 shown in Fig. 9共a兲. The data are constructed by averaging together 100 configurations. The curve shows fluctuations that decay with distance from the wall, eventually reaching a plateau. Using Eq. 共2兲, we obtain decay lengths ␭3D = 0.77ds and 0.73ds for the small and large particle curves, respectively. These length scales are similar to the length scales obtained in the 2D case 共␭2D = 0.85ds and 0.72ds for small and large particles兲. To compare all 3D ␳共y兲 distributions for different h, we construct the image representation used to compare 2D configurations in Fig. 5. The data for the 3D configurations are shown in Fig. 10. Again there are dark vertical strips arising from particles forming layers near the wall. Like in 2D, the density approaches the “bulk” rcp value far from the wall.

FIG. 10. An image representation comparing the number density distributions of 3D configurations for many different h. Black pixels represent a relative value of 1 and white represents a relative value of 0. A grayscale is used to represent relative values between 0 and 1. The pixel widths are 0.1ds horizontally and 0.2 vertically.

1 兺 Y 6m共␪ij, ␾ij兲. n jK j

共4兲

In the above equation m = 兵−6 , . . . , 0 , . . . , 6其, and thus qˆi,6 is a 13 element complex vector, which is assigned to every particle i in the system. The sum in Eq. 共4兲 is taken over the j nearest neighbors of the ith particle, n j is the total number of neighbors, and K is a normalization constant so that qˆi,6 · qˆi,6 = 1. For two particles i and j that are nearest neighbors, Y 6m共␪ij , ␾ij兲 is the spherical harmonic associated with the vector pointing from particle i to particle j, using the angles ␪ij and ␾ij of this vector relative to a fixed axis. Next, any two particles m and n are considered “ordered neigh* ⬎ 0.5 关80,81兴. Finally, we quantify the local bors” if qˆm,6 · qˆn,6 order within the system by the number of ordered neighbors Nb a particle has. Figure 9共b兲 is a plot of the average number of ordered neighbors particles have 具Nb典 as a function of distance y from the wall. In comparison with Fig. 9共a兲, this plot shows that local order is mostly seen within layers not between layers. Also we see that 具Nb典 converges to an asymptotic value of ⬇1.3, confirming that the system is disordered. 共Values of Nb ⱖ 8 are considered crystalline 关81兴.兲 We use Eq. 共2兲 to characterize a length scale for the decay in 具Nb典, giving ␭ = 1.3ds. The asymptotic limit of 具Nb典共y兲 in Fig. 9共b兲 agrees with the average value of Nb found for 15 large simulations with 2500 particles and periodic boundary conditions, confirming that the local structure in the confined case converges to the bulk rcp state far from the walls. Our results show that in both 2D and 3D, confinement induces changes in structural quantities near the walls, with a decay toward the “bulk” values characterized by length scales no larger than dl. The only prior work we are aware of with related results are a computational study 关40兴 and an experimental study 关41兴 of collections of monodisperse particles confined in a large silo. The simulation by Landry et al. primarily focused on the force network within the silo. They show one plot of the local packing fraction as a function of distance from the silo wall. Similar to our results, this local packing fraction showed fluctuations that decayed monotonically. In their paper, they state a decay length of ⬇4dl; however, it appears that they drew this conclusion by estimating the value by eye. Applying Eq. 共2兲 to their data, we find ␭ on the order of dl, close to the value found in our simulations. The experimental study by Seidler et al. reported on the local bond orientational order parameter, which showed oscillation that decayed with distance from the wall. They reported a decay length of ␭ ⬇ dl using an exponential fit. The length scales from these two studies are slightly larger than those found in our work. IV. MODEL

Our results for ␾共h兲 can be understood with a simple model incorporating an effective boundary layer and a bulk-

051305-7

PHYSICAL REVIEW E 80, 051305 共2009兲

KENNETH W. DESMOND AND ERIC R. WEEKS

FIG. 11. 共Color online兲 Illustration of the model for ␾共h兲. The model breaks a configuration with confining width h into three regions. The boundary layers are approximated to have a packing fraction ␾l and persist a distance ␦L into the sample, and the middle bulk region is approximated to have a packing fraction ␾rcp. These three parameters are assumed to be h independent.

like region. This model is an extension of one proposed by Verman and Banerjee 关31兴 and Brown and Hawsley in 1946 关32兴. In Fig. 11, we show a configuration of particles confined between two plates and divided into two boundary layers and a bulk region. The model of Refs. 关31,32兴 approximates the effect of the walls by assuming a lower effective volume fraction ␾l to the boundary layers. The central region is assumed to have a volume fraction ␾rcp, equal to the volume fraction for an infinite system. Of course, this model is an oversimplification that coarse grains the density near the walls, which in reality varies smoothly and nonmonotonically in space, as Sec. III demonstrates. Furthermore, this model will not capture the nonmonotonic behavior of Figs. 3 and 8, but it should capture the overall trend with h. In the original model of Refs. 关31,32兴, it was conjectured that the thickness of the boundary layer is ␦L = 1d for monotonic particles of diameter d. The experimental data they tested the model with were too limited to carefully check this assumption; here, we extend their model by allowing ␦L to be a free parameter. 共Clearly, our results, such as Fig. 7, confirm that ␦L ⬇ 1ds is a reasonable order of magnitude.兲 Using this simple model, ␾ can be approximated by the weighted average ␾ = h−2h␦L ␾rcp + 2␦hL ␾l 共in either 2D or 3D, with different values of the parameters depending on the dimension兲. Reducing this equation further, we obtain

␾ = ␾rcp −

C , h

FIG. 12. 共Color online兲 The upper black curve is a plot of ␾共1 / h兲 for 2D configurations and the red 共dark gray兲 line going through the curve is a fit from the model 关Eq. 共5兲兴. Likewise, the lower black curve is a plot of ␾共1 / h兲 for 3D configurations with the red 共dark gray兲 line going through the curve being another fit from the model.

overestimation on ␾rcp. When the data for both curves are fitted for h ⱖ 8 共1 / h ⬍ 0.125兲, the actual values for ␾rcp are obtained. The dipping of the ␾共h兲 curve below the line for large 1 / h is perhaps due to the layering each wall produces affecting the layering produced by the opposite wall 共see Figs. 5 and 10兲. Another possibility is that this reflects the breakdown of the model when h ⬇ 2␦L. That is, when the thickness of the sample is such that the two boundary layers begin to overlap, the model would not be expected to work. To provide further credence to the model, we also perform 2D rcp simulations with a fixed circular boundary or a fixed square boundary. Figure 13 shows a plot of ␾共h兲 for both the circular boundary 共green points兲 and the square boundary 共blue points兲. In analogy with our prior results, h is the wallto-wall distance: for the circular boundary, h is the diameter normalized by ds, and for the square boundary h is the side length L normalized by ds. As before with two parallel flat boundaries, we see that ␾ increases to an asymptotic limit.

共5兲

where we define the boundary packing parameter C = 2␦L共␾rcp − ␾l兲, which quantifies how the wall influences the packing fraction near the boundary. Note that this is the same form for ␾共h兲 obtained from considering a first-order correction in terms of 1 / h and is the same empirical form assumed by Scott 关5兴. We investigate the merit of this model by fitting the data to Eq. 共5兲, which only contains two fitting parameters ␾rcp and C. The data in both Figs. 3 and 8 are fitted to Eq. 共5兲. The fits are shown as the red lines in these earlier figures and also in Fig. 12, where the data are plotted as functions of 1 / h to better illustrate the success of this model. The fits give for 2D ␾rcp = 0.844 and C = 0.317 and for 3D ␾rcp = 0.646 and C = 0.233. Both fits give values for ␾rcp that are slightly larger, but not by much, than ␾rcp reported earlier in the paper. In Fig. 12, it can be seen that the packing fraction for large 1 / h dips significantly below the fitting line, due to the fluctuations in ␾共h兲 for small h; this is responsible for the

FIG. 13. 共Color online兲 ␾共h兲 for disks confined within a circular boundary or square boundary, as indicated. For the circular boundary, h is defined as the system’s diameter normalized by ds, and for the square boundary system h is defined as the system’s side length normalized by ds. The red 共light gray兲 curves are fits to the data using our model. The image at the lower right is an rcp configuration confined in a circular boundary with h = 21. Small particles are rendered as green 共medium gray兲 and large particles are rendered as blue 共dark gray兲.

051305-8

PHYSICAL REVIEW E 80, 051305 共2009兲

RANDOM CLOSE PACKING OF DISKS AND SPHERES IN…

␦L is a third parameter for the fit; fortunately, extending the model to a 3D case, where all directions are confined, would predict a cubic fit but without introducing a fourth parameter. V. CONCLUSION

FIG. 14. 共Color online兲 ␾共1 / h兲 for disks confined within a circular boundary or square boundary, as indicated; the data are the same as Fig. 13 along with the definitions of h. The red 共light gray兲 curves are fits to the data using our model. The black dashed lines are linear fits to the data. The image at the upper right is an rcp configuration confined in a square boundary with h = 15.5. Small particles are rendered as green 共medium gray兲 and large particles are rendered as blue 共dark gray兲.

Note that the data are noisier for two reasons. First, given our algorithm 共Sec. II, for samples that are confined in all directions兲, we can only choose the number of particles we start with; we have no control over the final system size when the particles jam. Due to random fluctuations, we can run the simulations many times with the same number of particles and each time find a different final value for h 共and ␾rcp兲. This limits our ability to sample enough data at a particular h to reduce the noise and/or look for nonmonotonic effects. Second, there are many fewer particles in these simulations, thus, reducing the statistics. Normally, this could be compensated by increasing the number of simulation runs, but the first problem 共lack of precise control over h兲 frustrates this. Adapting the model to a circular boundary with diameter h or to a square boundary with side length h, we find that both situations give

␾共h兲 = ␾rcp −

2C 2C␦L + , h h2

共6兲

where C = 2␦L共␾rcp − ␾l兲 as before. The data in Fig. 13 are fitted to this model, shown as the red lines. For the circular boundary, the fit gives ␾rcp = 0.846, C = 0.371, and ␦L = 1.51, and for the square boundary the fit gives ␾rcp = 0.848, C = 0.340, ␦L = 1.14. These fits give ␾rcp values close to the ␾rcp values reported earlier in the paper and C values similar, but slightly different, than that found for one fixed flat boundary. Interestingly, the fits give values of ␦L commensurate to the ␭ values previously computed, demonstrating that the boundary produces a thin boundary layer of about 1–2 characteristic particle sizes thick that is primarily responsible for lowering the global packing fraction. Finally, to demonstrate the quality of the fits we show a plot of ␾共1 / h兲 in Fig. 14. In this figure, the red line is the fit from Eq. 共6兲, while the black dashed line is a linear fit in 1 / h. Both fits are reasonable, and the data are not strong enough to determine which is better. We thus note only that the model suggests we should use the quadratic fit for these cases and that the values of ␦L obtained are reasonable ones.

In this paper, we have shown how a confining boundary alters the structure of random close packing by investigating simulated rcp configurations confined between rigid walls in 2D and 3D. We find that confinement lowers the packing fraction and induces heterogeneity in particle density, where particles layer in bands near the wall. The structure of the local packing decays from a more ordered packing near the wall to a less ordered packing in the bulk. All measures of local order and local density decay rapidly to their bulk values with characteristic length scales on the order of particle diameters. Thus, the influence of the walls is rapidly forgotten in the interior of the sample, with confinement having the most notable effects when the confining dimension is quite small, perhaps less than 10 particle diameters across. The results are well fit by a three-parameter model dating back to 1946 关31,32兴, with our results suggesting that the third parameter 共an effective boundary layer thickness兲 should be a free parameter rather than constrained. To first order, this model suggests that the primary influence of the boundary is quantified by one parameter C, which is the product of a length scale and a volume fraction reduction. This parameter, the boundary packing parameter, thus, quantifies the overall influence of a boundary, near that boundary. Since the model assumes nothing about shape, this model should equally apply to other geometries as well. These findings have implications for experiments investigating the dynamics of densely packed confined systems 共i.e., colloidal suspensions or granular materials兲. For example, our work shows that for small h the packing fraction has significant variations at small h 共most clearly seen in 2D, for example, Fig. 3兲. For dense particulate suspensions with ␾ ⬍ ␾rcp, flow is already difficult. By choosing a value of h with a local maximum in ␾rcp共h兲, a suspension may be better able to flow, as there will be more free volume available. Likewise, a poor choice of h may lead to poor packing and enhanced clogging. A microfluidic system with a tunable size h may be able to vary the flow properties significantly with small changes in h, but our work implies that control over h needs to be fairly careful to observe these effects. Of course, these effects will be obscured by polydispersity in many systems of practical interest; however, our work certainly has implications for microfluidic flows of these sorts of materials, once the minimum length scales approach the mean particle size. Our work has additional implications for experiments on confined glasses 关50,51,65,67–72兴. As mentioned in the introduction, confinement changes the properties of glassy samples, but it is unclear if this is due to finite-size effects or due to interfacial influences from the confining boundaries 关73兴. Our results show that dense packings have significant structural changes near the flat walls, suggesting that indeed interfacial influences on materials can be quite strong at very short distances, assuming that the structural changes couple

051305-9

PHYSICAL REVIEW E 80, 051305 共2009兲

KENNETH W. DESMOND AND ERIC R. WEEKS

ACKNOWLEDGMENTS

with dynamical behavior. Furthermore, the nonmonotonic behavior of ␾rcp that we see suggests that experiments studying confined glassy materials could see interesting nonmonotonic effects, if the sample thickness can be carefully controlled.

We thank C. S. O’Hern for helpful conversations. This work was supported by the National Science Foundation under Grant No. DMR-0804174.

关1兴 A. E. R. Westman and H. R. Hugill, J. Am. Ceram. Soc. 13, 767 共1930兲. 关2兴 H. E. White and S. F. Walton, J. Am. Ceram. Soc. 20, 155 共1937兲. 关3兴 A. H. Boerdijk, Philips Res. Rep. 7, 303 共1952兲. 关4兴 D. P. Haughey and G. S. G. Beveridge, Can. J. Chem. Eng. 47, 130 共1969兲. 关5兴 D. G. Scott, Nature 共London兲 188, 908 共1960兲. 关6兴 R. M. German, Particle Packing Characteristics 共Metal Powder Industries Federation, Princeton, NJ, 1989兲. 关7兴 S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 共2000兲. 关8兴 J. D. Bernal and J. Mason, Nature 共London兲 188, 910 共1960兲. 关9兴 O. K. Rice, J. Chem. Phys. 12, 1 共1944兲. 关10兴 W. O. Smith, P. D. Foote, and P. F. Busang, Phys. Rev. 34, 1271 共1929兲. 关11兴 S. F. Edwards, Granular Matter 共Springer-Verlag, Berlin, 1994兲. 关12兴 C. Radin, J. Stat. Phys. 131, 567 共2008兲. 关13兴 M. Jerkins, M. Schröter, H. L. Swinney, T. J. Senden, M. Saadatfar, and T. Aste, Phys. Rev. Lett. 101, 018301 共2008兲. 关14兴 R. Pal, Polym. Eng. Sci. 48, 1250 共2008兲. 关15兴 G. Lois, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev. Lett. 102, 015702 共2009兲. 关16兴 R. Zallen, Physics of Amorphous Solids 共Wiley, New York, 1983兲. 关17兴 C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 70, 043302 共2004兲. 关18兴 R. K. Mcgeary, J. Am. Ceram. Soc. 44, 513 共1961兲. 关19兴 A. R. Kansal, S. Torquato, and F. H. Stillinger, J. Chem. Phys. 117, 8212 共2002兲. 关20兴 R. Al-Raoush and M. Alsaleh, Powder Technol. 176, 47 共2007兲. 关21兴 K. Lochmann, L. Oger, and D. Stoyan, Solid State Sci. 8, 1397 共2006兲. 关22兴 H. J. H. Brouwers, Phys. Rev. E 74, 031309 共2006兲. 关23兴 T. Okubo and T. Odagaki, J. Phys.: Condens. Matter 16, 6651 共2004兲. 关24兴 P. Richard, L. Oger, J. P. Troadec, and A. Gervois, Eur. Phys. J. E 6, 295 共2001兲. 关25兴 R. L. Brown and P. G. W. Hawksley, Nature 共London兲 156, 421 共1945兲. 关26兴 H. S. Leftwich, Nature 共London兲 156, 753 共1945兲. 关27兴 A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato, and P. M. Chaikin, Science 303, 990 共2004兲. 关28兴 K. Desmond and S. V. Franklin, Phys. Rev. E 73, 031306 共2006兲. 关29兴 Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 100, 245504 共2008兲.

关30兴 P. C. Carman, Trans. Inst. Chem. Eng. 15, 150 共1937兲. 关31兴 L. C. Verman and S. Banerjee, Nature 157, 584 共1946兲. 关32兴 R. L. Brown and P. G. W. Hawksley, Nature 共London兲 157, 585 共1946兲. 关33兴 G. D. Scott and D. M. Kilgour, J. Phys. D 2, 863 共1969兲. 关34兴 K. Ridgway and K. J. Tarbuck, Br. Chem. Eng. 12, 384 共1967兲. 关35兴 B. R. Aïm and L. P. Goff, Powder Technol. 1, 281 共1968兲. 关36兴 A. G. Dixon, Can. J. Chem. Eng. 66, 705 共1988兲. 关37兴 M. Suzuki, T. Shinmura, K. Iimura, and M. Hirota, Powder Technol. 19, 183 共2008兲. 关38兴 R. Zou, Chem. Eng. Sci. 50, 1504 共1995兲. 关39兴 R. P. Zou and A. B. Yu, Chem. Eng. Sci. 51, 1177 共1996兲. 关40兴 J. W. Landry, G. S. Grest, L. E. Silbert, and S. J. Plimpton, Phys. Rev. E 67, 041303 共2003兲. 关41兴 G. T. Seidler, G. Martinez, L. H. Seeley, K. H. Kim, E. A. Behne, S. Zaranek, B. D. Chapman, S. M. Heald, and D. L. Brewe, Phys. Rev. E 62, 8175 共2000兲. 关42兴 L. Vanel, P. Claudin, J. P. Bouchaud, M. E. Cates, E. Clément, and J. P. Wittmer, Phys. Rev. Lett. 84, 1439 共2000兲. 关43兴 U. Marconi, Physica A 280, 279 共2000兲. 关44兴 P. Claudin and J.-P. Bouchaud, Phys. Rev. Lett. 78, 231 共1997兲. 关45兴 J. Landry, Powder Technol. 139, 233 共2004兲. 关46兴 D. M. Mueth, Phys. Rev. E 67, 011304 共2003兲. 关47兴 J. Ritvanen and P. Jalali, Physica A 387, 5381 共2008兲. 关48兴 L.-W. Teng, P.-S. Tu, and L. I, Phys. Rev. Lett. 90, 245004 共2003兲. 关49兴 C. Murray, MRS Bull. 23, 33 共1998兲. 关50兴 Z. T. Németh and H. Löwen, Phys. Rev. E 59, 6824 共1999兲. 关51兴 M. Alcoutlabi and G. B. McKenna, J. Phys.: Condens. Matter 17, R461 共2005兲. 关52兴 P. Gallo, M. Rovere, and E. Spohr, J. Chem. Phys. 113, 11324 共2000兲. 关53兴 S. Granick, Science 253, 1374 共1991兲. 关54兴 J. R. Henderson, Mol. Phys. 105, 2345 共2007兲. 关55兴 J. Mittal, V. K. Shen, J. R. Errington, and T. M. Truskett, J. Chem. Phys. 127, 154513 共2007兲. 关56兴 J. Mittal, J. R. Errington, and T. M. Truskett, J. Chem. Phys. 126, 244708 共2007兲. 关57兴 J. Mittal, T. M. Truskett, J. R. Errington, and G. Hummer, Phys. Rev. Lett. 100, 145901 共2008兲. 关58兴 I. Cohen, T. G. Mason, and D. A. Weitz, Phys. Rev. Lett. 93, 046001 共2004兲. 关59兴 K. To, P.-Y. Lai, and H. K. Pak, Physica A 315, 174 共2002兲. 关60兴 I. Zuriguel, L. A. Pugnaloni, A. Garcimartín, and D. Maza, Phys. Rev. E 68, 030301共R兲 共2003兲. 关61兴 S. Redner and S. Datta, Phys. Rev. Lett. 84, 6018 共2000兲. 关62兴 S. B. Fuller, E. J. Wilhelm, and J. M. Jacobson, J. Microelectromech. Syst. 11, 54 共2002兲.

051305-10

PHYSICAL REVIEW E 80, 051305 共2009兲

RANDOM CLOSE PACKING OF DISKS AND SPHERES IN… 关63兴 K. Sharp and R. Adrian, Microfluid. Nanofluid. 1, 376 共2005兲. 关64兴 G. M. Whitesides, Nature 共London兲 442, 368 共2006兲. 关65兴 C. R. Nugent, K. V. Edmond, H. N. Patel, and E. R. Weeks, Phys. Rev. Lett. 99, 025702 共2007兲. 关66兴 P. N. Pusey and W. van Megen, Nature 共London兲 320, 340 共1986兲. 关67兴 D. Morineau, Y. Xia, and C. A. Simionesco, J. Chem. Phys. 117, 8966 共2002兲. 关68兴 J. Schüller, Y. B. Mel’nichenko, R. Richert, and E. W. Fischer, Phys. Rev. Lett. 73, 2224 共1994兲. 关69兴 C. L. Jackson and G. B. McKenna, J. Non-Cryst. Solids 131133, 221 共1991兲. 关70兴 K. Kim and R. Yamamoto, Phys. Rev. E 61, R41 共2000兲. 关71兴 P. A. Thompson, G. S. Grest, and M. O. Robbins, Phys. Rev. Lett. 68, 3448 共1992兲. 关72兴 P. Scheidler, W. Kob, K. Binder, and G. Parisi, Philos. Mag. A 82, 283 共2002兲.

关73兴 F. He, L. M. Wang, and R. Richert, Eur. Phys. J. Spec. Top. 141, 3 共2007兲. 关74兴 N. Xu, J. Blawzdziewicz, and C. S. O’Hern, Phys. Rev. E 71, 061306 共2005兲. 关75兴 A. S. Clarke and J. D. Wiley, Phys. Rev. B 35, 7350 共1987兲. 关76兴 M. Matsumoto and T. Nishimura, ACM Trans. Model. Comput. Simul. 8, 3 共1998兲. 关77兴 J. Nocedal and S. J. Wright, Numerical Optimization 共Springer, New York, 1999兲. 关78兴 A. Donev, S. Torquato, F. H. Stillinger, and R. Connelly, Phys. Rev. E 70, 043301 共2004兲. 关79兴 A. H. Marcus and S. A. Rice, Phys. Rev. Lett. 77, 2577 共1996兲. 关80兴 P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 共1983兲. 关81兴 U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey, and D. A. Weitz, Science 292, 258 共2001兲.

051305-11