Three-Dimensional Monte Carlo Simulations of Density Relaxation †

The compaction and densification of granular matter is an issue of importance in the industrial sector in that an increase in packing efficiency in the processing of bulk solids is often desirable in reducing costs and in meeting consumer demands. However, in some situations, the compaction resulting from tapping or vibrations in handling operations and long term transport may be detrimental to the quality of the end product. Beyond these practical matters, the ability of a bulk solid to undergo a change in its bulk density due to external disturbances is a fundamental attribute of granular materials that is not well-understood. Indeed, the factors that have a bearing on the phenomenon (e.g., particle properties, containment geometry, environmental conditions, and the nature of the imposed disturbances) are rather broad and consequently the relationship between them requires further study. The increase in bulk density experienced by contained granular materials subjected to taps and/or continuous vibrations (often referred to as ‘density relaxation’) has been well-documented in the early literature 6, 17-19, 33, 49, . From a broader perspective, research on this phenomenon has its foundations in the many fundamental studies on the packing of particles (for example, see 1, 3, 9, 10, 14, 16, 21, 25, 26, 40, 46-48, 53, 54)). In recognition of the importance of developing a better understanding of the nature of the granular state, and motivated in part by recent experiments 27, , density relaxation has received renewed interest as evidenced by the recent upsurge in both modeling and computational studies, such as 2, 11, 13, 20, 24, 29-32, 34, 42, . In this paper, the effect of discrete taps applied to a vessel of granular material is modeled using a Monte Carlo approach in three-dimensions. This methThree-Dimensional Monte Carlo Simulations of Density Relaxation


Introduction
The compaction and densification of granular matter is an issue of importance in the industrial sector in that an increase in packing efficiency in the processing of bulk solids is often desirable in reducing costs and in meeting consumer demands.However, in some situations, the compaction resulting from tapping or vibrations in handling operations and long term transport may be detrimental to the quality of the end product.Beyond these practical matters, the ability of a bulk solid to undergo a change in its bulk density due to external disturbances is a fundamental attribute of granular materials that is not well-understood.Indeed, the factors that have a bearing on the phenomenon (e.g., particle properties, containment geometry, environmental conditions, and the nature of the imposed disturbances) are rather broad and consequently the relationship between them requires further study.
The increase in bulk density experienced by contained granular materials subjected to taps and/or continuous vibrations (often referred to as 'density relaxation') has been well-documented in the early literature 6, 17-19, 33, 49, 50) .From a broader perspective, research on this phenomenon has its foundations in the many fundamental studies on the packing of particles (for example, see 1, 3, 9, 10, 14, 16, 21, 25, 26,  40, 46-48, 53, 54)).In recognition of the importance of developing a better understanding of the nature of the granular state, and motivated in part by recent experiments 27,[37][38][39] , density relaxation has received renewed interest as evidenced by the recent upsurge in both modeling and computational studies, such as 2, 11, 13, 20, 24, 29-32, 34, 42, 52) .
In this paper, the effect of discrete taps applied to a vessel of granular material is modeled using a Monte Carlo approach in three-dimensions.This meth-

Three-Dimensional Monte Carlo Simulations of Density Relaxation †
Oleksandr Dybenko and Anthony D. Rosato* Granular Science Laboratory, Dept. of Mechanical Engineering, New Jersey Institute of Technology odology does not consider the detailed dynamical interactions between particles, as would be the case with dissipative molecular dynamics simulations 56) .Rather, the effect of the intensity of the taps on the equilibrium bulk solids fraction is examined.In this regard, the intensities that were selected roughly correspond to relatively large accelerations of the vessel.A broad set of tap intensities was used, thereby allowing us to obtain a correspondingly wide range of solids fractions, from a loose configuration, rather like a 'poured' assembly, to a relatively dense structure with local crystalline order.The vessel itself was chosen to have a smooth (i.e., not bumpy) floor so as not to preclude the well-known 'ordering' effect of a flat plane on packing of uniform spheres 23,41,43) .The outline for the remainder of this paper is as follows.The next section presents the details of the simulation technique that was developed to evolve the system from an initial low solids fraction state to a denser configuration.We describe the means by which taps are applied to the system of particles, the physical interpretation of these taps, the computation of equilibrium averages and the convergence criterion used.Simulation results are presented in Section 3. Here, the progression of the bulk solids fraction with tap number is shown for various tap intensities, in addition to a study of the equilibrium solids fraction as a function of the tap intensity.These simulation results are also compared with various models from the literature.In Section 4, results on the microstructure are given and quantified through the contact number distribution and radial distribution function.Finally, the conclusions are presented in Section 5.

Simulation Methodology
In this section, the approach used to generate equilibrium assemblies of spheres is described.The quantity of interest is the bulk solids fraction ν, which is calculated from the simulation data as the ratio of the total sphere volume to the volume occupied.In carrying out the calculations of ν, the upper 20% of the packing was not included so as to minimize the influence of the top surface and to simplify the numerical procedure without a loss in accuracy.We remark that the heap formed 8,15) at the surface as a consequence of a granular material's angle of repose in a physical experiment was not a phenomenon that could be captured by our simulation method.
A Monte Carlo scheme was used to simulate the density relaxation process that occurs as a result of a series of discrete taps applied to a 'box' containing N spherical particles of diameter d.The computational domain (i.e. the 'box') is a rectangular parallelepiped having periodic boundar y conditions in the lateral directions, which are hereby denoted as x and z.The lateral dimensions of the box (12d×12d) are chosen to optimize CPU time.In order to model the effect of the taps on the particles, each sphere is subjected to a height-dependent, vertical displacement whose magnitude is controlled by the lift intensityγ.This expansion of the system allows a greater portion of the bulk mass to reside near the floor of the box.For each selected value ofγ, the bulk solids fraction is traced as a function of the number of taps, and its equilibrium value ν ∞ (γ) is computed.The next section presents a detailed description of the Monte Carlo method.

Generation of the poured assembly
The simulation is initiated by randomly placing within the box N spheres that are subsequently allowed to collapse vertically under gravity to generate a configuration corresponding to a loose ('poured') assembly.The spheres are 'hard'; that is, in order to avoid overlap of the spheres, for each i≠j, the distance | → ri− → rj| between the centers → ri =(xi, yi, zi) and rj =(xj, yj, zj) of particles is always greater than d+ε, where ε=10 −4 d .The Metropolis 35) algorithm is then used to advance the system to a stable configuration.A single Monte Carlo (MC) step in the Metropolis algorithm consists of the random selection of a particle (xi, yi, zi), followed by its assignment to a trial position in accordance to the prescription, Here, ξx, ξy, ξz, are random numbers sampled from a uniform distribution on [0, 1], and δ is a maximum allowed displacement.The procedure for selecting δ is described in the next paragraph.The trial position is accepted unconditionally as the new center coordinate if the change in the system potential energy is negative, i.e., ∆E ≡ mg Otherwise, (ΔE≥0) the trial position is accepted with probability e −βΔE .For the macroscopic particles that we are studying, βis very large 45) so that for all practical purposes, only downward displacements are allowed.Another particle is then selected at random, and the above procedure is repeated.
As the process described above advances through many thousands of MC steps, the space between adjacent particles is reduced and the bulk density increases.If the maximum allowed displacementδ were set to a small fixed value δ/d<<1 at the outset, the evolution of the system to a local equilibrium would be prohibitively slow.The situation would be the same if δ/d-1, since the percentage of accepted steps would be drastically reduced as the system density increases.In order to prevent this from taking place, δis dynamically changed from a maximum valueδ=d at the beginning of the process to a minimum δ=10 −3 d, with its reduction determined by the accepted fraction of steps every 10 4 MC steps.If more than half of the MC steps have been accepted, δremains unchanged; other wiseδis modified as follows:δ′=0.995δ.While other algorithms in the literature may be more computationally efficient 12,51) , this approach has a great deal of virtue in its simplicity.
The change in the bulk solids fraction is computed every one million MC steps and the Metropolis algorithm is terminated when where ν(n;γ) is the bulk solids fraction at tap number n and lift intensityγ.To ensure that this stopping criterion did not cause premature termination, in many test runs we required the criterion in (2) to be satisfied twice before ending the process.In none of these instances was a premature termination observed, allowing us to satisfy (2) only once for many of the results presented below.In all of the case studies reported in this paper, the initial random assembly collapsed to a loose configuration having a bulk solids fraction in the range 0.56 to 0.58.In general, the poured assemblies filled the box to a height H/d 11 for the system size considered in this work.

Tapping procedure
The poured assembly is dislodged from its local energy minimum by applying vertical position dependent displacements to each particle in accordance to The factor γ represents the effect of a single tap so that in a qualitative sense, a greater impulse applied to the container corresponds to a larger value ofγ.We remark that the form (3) for vertical displacements is equivalent to that used by Barker and Mehta 7) , who also imposed random lateral displacements on the particles at each tap.The influence of applying this additional lateral displacement on the equilibrium solids fraction is considered later in this paper.After expanding the vertical coordinates using (3), the system is allowed to settle using the Metropolis algorithm as described in the previous section.Satisfaction of the stopping criterion (2) to signal completion of a single 'tap' applied to the system typically required 17 to 20 million MC steps (~ 40 CPU seconds on our computers).Fig. 1 shows a sequence of configurations, starting from the initial random placement, the 'poured assembly', its expansion via a single tap (3), and its state in accordance with the criterion (2).Note that the shading in this figure is only for the purpose of visualization.

Equilibrium bulk solids fraction
A tapping cycle is initiated with the poured particle assembly as described in Section 2.1.A single experiment (or realization) consists of a sequence of individual taps (Section 2.2).As used here and henceforth, the term 'tap' refers to the expansion by (3) and subsequent relaxation.Each realization is carried out using a fixed lift intensityγ; the total number of taps to obtain equilibrium depends onγ.For smaller values ofγ, more taps were needed to reach an equilibrium solids fraction.Typical plots of bulk solids fraction against the number of taps NT are shown in Poured assembly (ν＝0.568),(c).Configuration after applying one tap, and (d).Final configuration at 1000 taps.Note that particles are shaded only for the purpose of visualization.
Fig. 2 for two different realizations atγ=1.25.The graphs exhibit an overall increase in the solids fraction with the number of taps, a general trend that was evident in all realizations.The convergence to the equilibrium configuration also exhibited dependence upon the choice ofγ, which will be demonstrated in subsequent figures.
The ensemble averaged bulk solids fraction 〈ν(n;γ)〉 was computed over M realizations.In order to identify the minimum number of realizations M so that a statistically correct equilibrium value was obtained, the standard deviation 〈σ(n;γ)〉 as a function of the tap number n was determined as where The behavior of 〈σ(n;γ)〉 was monitored by increasing the number of realizations M until, over a sufficient number of taps NT, 〈σ(NT;γ)〉< 0.001 (6)   which essentially defines equilibrium in the simulations.The behavior of the standard deviation 〈σ(NT;1.25)〉depicted in Fig. 3 demonstrates that criterion ( 6) is reached when of NT-800 for this value ofγ.The ensemble averaged solids fraction as a function of M is shown in Fig. 4 for the same case when NT=1000.The vertical lines in the figure denote error bars in the ensemble averaged solids fraction.In general, it was found that M 100 realizations were adequate to satisfy the criterion in (6).

Simulation Results
The value of the equilibrium bulk solids fraction depends on the lift intensity γ, which is roughly analogous to changing the intensity of the taps on the containment vessel in a physical experiment.As will be seen, the results demonstrate that it is possible to obtain a rather wide range of equilibrium solids fractions, from loosely packed assemblies to rather  dense configurations whose radial distribution functions contains peaks associated with the presence of crystalline structures.It is noted however that these dense configurations are largely promoted by the ordering influence of the plane floor -a phenomenon that we did not wish to preclude by using a bumpy base.
For the pouring phase of the simulation, we obtained 'poured' solids fractions in the range 0.560 to 0.580, corresponding to a random loose structures reported in the literature 27,37,40,44) .Starting from the poured structure, taps are applied to the system following the procedure described in Sec.2.2 and Sec.2.3.The ensemble-averaged bulk solids fraction 〈ν(n;γ)〉 for lift intensities γ between 1.25 and 1.50 as a function of the tap number is shown in Fig. 5.
For the purpose of clarity, only every twentieth data point is plotted; here each symbol represents the average bulk solids fraction over 100 realizations.It can be seen in the figure that taps of smaller intensity eventually yield larger equilibrium bulk densities.Simulations were also performed at γ＝1.2; however these results were not included in Fig. 5 because of the larger number of taps (NT ≥2,000) required to reach equilibrium.In fact, for γ＝1.15, we carried out 10,000 taps, which were not sufficient to satisfy the equilibrium criterion in (6).This is not unexpected given the trend shown in Fig. 6 for the number of taps required to reach 99.5% of the equilibrium solids fraction.The deviations (equations ( 4) and ( 5)) for γ ＝1.25, 1.35 and 1.50 as a function of the tap number in Fig. 7 exhibit a clear-cut dependence on lift intensity γ.The peak observed in the γ＝1.25 case was observed to be a characteristic feature at small values of the intensity.However, as the lift intensity became smaller, the peak in the deviation grew so that it re-quired significantly more taps to satisfy criterion (6).At the smaller intensities γ considered here, the system tends to visit a greater portion of the phase space in approaching its final lower energy state.
In order to determine if random lateral displacements of the particles in sync with vertical taps would significantly change the equilibrium bulk solids fraction, we modified our Monte Carlo procedure by incorporating this feature of Mehta et al.'s algorithm 7) .That is, in conjunction with the vertical motion prescribed by (3), each particle is assigned a random lateral displacement within a disk of radius ρ given by, The parameter ρ is computed at each tap as the average (over Np particles) of the minimum nearneighbor lateral distance, Here, , where Fig. 7 The standard deviation of bulk solids fraction〈σ(n;γ)〉versus tap number N T for tap intensity γ＝1.25 (Δ), 1.35 (♢) and 1.50 (*).Smaller values of γ yield higher variances.
Fig. 6 The number of taps required to attain 99.5% of equilibrium bulk solids fraction versus the tap intensity γ.A reduction of γ results in a greater number of taps to reach equilibrium.the minimum is taken over all par ticles j whose centers are within two diameters from particle k.In Fig. 8, we compare the bulk solids fraction versus tap number at lift intensity γ＝1.25 with and without application of the random lateral displacements (7).
The trends of the two graphs are very similar, and we observed steady-state values of the solids fraction that were statistically equivalent.A further check was done by introducing lateral displacements in a system that had equilibrated without them in order to determine if it would return to the same equilibrium solids fraction.In the results of this study, not shown here for the sake of brevity, it was observed that random lateral displacements did not move the system out of equilibrium.

Evolution behavior
A phenomenological relationship for the evolution observed in the Monte Carlo simulations is described in this section.This relation proved to be quite useful in the prediction of the equilibrium packing fractions for the case studies presented in this paper.However, it must be noted that the connection between the Monte Carlo evolution and that which occurs in dynamic simulation is not clearly understood and is a topic that requires further investigation.Nonetheless, a comparison of the evolution of the bulk solids fraction computed from our simulations is made with relationships reported in the literature.
Our simulation data could be fit to a hyperbolic tangent model (solid lines in Fig. 5) of the form where γ 0 (γ) and γ ∞ (γ) are its initial and equilibrium values, respectively and A(γ) is the controlling rate constant.The use of ( 9) is motivated in part by the rate law dρ(t) ＝ − β ρ dt , the solution of which suggests a difference equation of the form, The general solution of this difference equation is given by (11)   By expressing the model in (9) in terms of Linz's 30) compaction ratio ρ n ≡ The solution (11) with B＝1, β＝2A C＝2(1− e 2A ) is actually the first two terms of the asymptotic expansion of (12).Parameters determined from the least-squares fit to (9) using the simulation results appear in Table 1 for γ≥1.2.Recall that these parameters are the initial and final solids fractions ν 0 (γ) and ν ∞ (γ), respectively, and the rate constant A(γ).Observe that the initial solids fractions ν 0 (γ) from the fit fall within the range 0.560 to 0.580 generated by the simulation.The graph of ν ∞ (γ) versus γ in Fig. 9 shows that a reduction in the lift intensity results in the creation of a more dense system, but at the expense of an increasingly growing number of taps (see Fig. 6) to reach equilibrium via criterion (6).Values of ν ∞ (γ) from the fit closely matched that obtained from the simulation for γ≥1.2.Furthermore, the model dis-  played rather good agreement with the progression of the data with tap number n, although (as expected) less than perfect conformity was found at the onset of the densification process.However, the hyperbolic tangent function was not observed fit the lower intensity data (at γ＝1.1 and 1.15) as well since computational limitations precluded the generation of a data set of sufficient size in the sense of being close to equilibrium.In fact, a crude backward extrapolation of the data in Fig. 6 suggested that for lift intensity γ＝1.1 upwards of 40,000 taps may be required.
For the sake of completeness, we have also considered fits of the simulation data to other models in the literature.In particular, comparisons were done with the stretched exponential 28,55) , and the inverse log law 27) .The stretched exponential takes the form There are four parameters in this model, corresponding to the initial solids fraction ν 0 , the equilibrium solids fraction ν ∞ , the Monte Carlo time scale τ and the stretching exponent β 28) .The fits to the data appear in Table 2 from those cases which equilibrium was attained in accordance with ( 6); these fits do not reveal a consistent value or a trend in the stretching exponent β.Furthermore, neither τ nor τ β revealed a decrease with intensity γ, as would be expected for a time scale measure.Despite these issues, the fits matched the data for the tap intensities considered in this study, with statistical measures of goodness of fit (such as R 2 values) not substantially different than those seen in the hyperbolic tangent fits.
The inverse log model is given by where the four parameters are the initial and equilibrium solids fractions ν 0 and ν ∞ , respectively, the Monte Carlo time scale τ and an amplification factor B. For tap intensities γ ≥1.2, poor correlation between this model and the simulation data was found since equilibrium was attained more rapidly than could be accommodated by the form (14).In the lower tap intensity regime (i.e.,γ＝1.1 and 1.15), Table 3 shows that nearly constant values of the ratio B/τ were obtained as a function of the number of taps, with the value of this ratio being dependent on γ.The possibility that B/τ is a relevant parameter for low tap intensities is also suggested by the experimental data in Fig. 5 of Knight et al. 27) .An exact correlation between γ in the Monte Carlo simulations and the bed's response to the sinusoidal taps, as in the experiments of Knight et al., is a matter of further research.However, rough scaling arguments (done by equating the normalized average bed expansion in the simulations with the experimental tap amplitude α) suggest that the lower intensities (γ＝1.1 and 1.15) correspond to accelerations Г≡αω 2 /g greater than the experimental values, but within the same order of magnitude.A similar fit analysis of the data to (14) for γ≥1.2 revealed that B/τ was not constant as a function of the number of taps; for example, the ratio increased monotonically by an order of magnitude for γ＝1.25 when 600 ≥ NT ≥ 3000.Further exploration of the significance of B/τ was carried out by fitting the low intensity data (γ＝1.1 and 1.15 ) to the functional form The results are given in Table 4, where it is observed that the values of D are close to B/τ obtained from the inverse log fits (Table 3).Furthermore, the increase of D with γ is in line with the expected behavior of a time scale parameter.These findings suggest the ratio B/τ is an important factor that captures the evolution of the data as a function of Monte Carlo time (i.e., tap number n).Note that (15)  is simply the first order approximation to the denominator of ( 14) with the identification that D≡B/τ.

Characterization of the Microstructure
The dependence of the configuration microstruc- ture was examined by computing ensemble-averaged (i.e., equilibrium) distribution of near-neighbor contacts as well as the radial distribution function.These computations were done for a variety of lift intensities.Distribution functions were determined by averaging on configurations in which the upper 20% of the particles was removed from the data so as to minimize the influence of the top surface and to simplify (without a loss in accuracy) the numerical procedure.
Spheres were defined as being in contact if the distance between their centers is less than or equal to 1.05d.The number of nearest neighbors is commonly referred to as the coordination number.Fig. 10 shows the expected result that the near-neighbor contact distributions shift to the right corresponding to an increase in mean coordination number〈NC〉with bulk solids fraction ν ∞ as determined from the fit to equation (9).A comparison of the simulated values 〈NC〉as a function of ν ∞ with the experimental measurements of Aste et al. 4) is shown in Fig. 11, where the simulation results are denoted by solid triangles.
Here, good agreement is observed with the trend of the data, as indicated by the linear regression fits appearing in the figure and a quantitative comparison of the data with the reported experimental averages.
The radial distribution function 22) is given by , where N is the number of particles, ρ is the bulk density (i.e., number of particles per unit volume where ν＝ρπd 3 /6 ), rij is the distance between the centers of particles i and j, and δ is the standard delta distribution.Thus, g − (r) may be considered as that quantity which, when multiplied by the bulk density, yields the local density.Computationally, g − (r) is approximated as an ensemble-average quantity taken at the last Monte Carlo time step (or tap).
where〈・〉represents the ensemble average, N(r) is the number of sphere centers within a spherical shell of radius r and thickness Δr, and ν is the bulk solids fraction.Due to limited computational facilities, more refined 'time' and ensemble averaging was not possible over the range of tap intensities considered in this study.However, several preliminary calculations indicated little difference in the radial distribution results when 'time' averaging was included.Fig. 12 shows〈g − (r/d)〉of the initial configuration and the result of three tap intensities corresponding to decreasing solids fractions ν ∞ .As the density increases, observe that the peaks in the distribution function become more pronounced.For ν ∞ ＝0.675, a power-law singularity at r＝d of the form ḡ(r) ≈ c 0 |r − r 0 | −α (17)   was found.The least-squares fit to (17) α＝ 0.528, in agreement with the value reported in the molecular dynamics studies of Silbert et al. 47) and later experiments of Aste 5) .The second and third peaks in the distribution appear at r/d ＝ √ 2 and√ 3. The influence of the plane floor in the simulations is enhanced by the taps, as can be seen in the radial distribution function of Fig. 12 at ν ∞ ＝0.675.At the lowest tapped solids fraction (circled in Fig. 12, for ν ∞ ＝0.583), it is observed that the two peaks at r/d -1.8 and 2.0 provide evidence of a local structural arrangement similar to the experimental results of Aste 5) .This structure does not appear in the poured assembly ν0＝0.567 as seen in the plot at the lower right hand corner of Fig. 12.We remark that this local arrangement differs from that seen by Nicolas et al. 36) , who obser ved hexagonal close-packed structures resulting from the application of horizontal cyclic shearing.
It is well-known that the presence of a flat plane boundary affects the local packing structure.In our simulations, this effect is enhanced through the taps as can be seen in Fig. 13, where the solids fraction as a function of the height is shown for the case ν ∞ ＝0.654.The graph was obtained by computing solids fraction within successively increasing layers starting from the floor.The oscillations are the signatures of a microstructure that is strongly affected by the presence of a planar floor.As the top of the packing is approached, there is the expected decay in the magnitude of the oscillations and eventual rapid decrease.Preliminary case studies of larger systems have demonstrated substantially the same effects.Further investigations of these phenomena are in progress and the results will be reported elsewhere.

Summar y and Conclusions
Density relaxation is a phenomenon of great interest because of its importance in the industrial section concerned with the handling and processing of bulk solids.In this paper, a Monte Carlo study was carried out to investigate the density relaxation process arising from discrete taps applied to a granular material consisting of monodisperse hard spheres within a laterally periodic computational volume.A single tap of intensity γ is modeled by applying vertical position-dependent displacements to the particles, which are then allowed to collapse under gravity via the Metropolis algorithm, thereby generating stable configurations.The algorithm is tailored so that motions which increase the system potential energy occur very rarely.The application of random lateral displacements in conjunction with the vertical expansion produced no effect on the final outcome.In order to accelerate the process and avoid trapping at local minima, the size of the neighborhood in which a particle is allowed to move was dynamically reduced as the system density increases.For tap intensity γ, the equilibrium value of the bulk solids fraction 〈ν(n;γ)〉 at tap n was computed as an ensemble average over a sufficient number of realizations to ensure a standard deviation of less than 0.001.It was found that more dense systems were created at smaller tap intensities, but this was at the expense of an increasingly growing number of taps to reach equilibrium.The evolution of the simulated data fit well to a hyperbolic tangent model of the form ν(n;γ) ＝ν0(γ) ＋ ［ν ∞ (γ) −ν0(γ)］tanh(A(γ)n) for tap intensities 1.2≤γ≤1.5.For these energetic taps, poor correlation between the inverse log model and the simulation data was because equilibrium was attained more rapidly than could be accommodated by this functional form.However, at smaller intensities that roughly corresponded to the magnitudes of the tap accelerations used in experiments 27) , the inverse log model was consistent with the data.The fit of the data to this model and to the form indicated that the ratio D -B/τ is a relevant parameter for low tap intensities.A fit to the stretched exponential ν(n)＝ν ∞ −(ν ∞ −ν0)e − (n/τ)β showed reasonably good agreement with the data, albeit neither τ nor τ β exhibited an expected decrease with intensity γ compatible with a time scale parameter.Quantification of the microstr ucture via the radial distribution function g(r/d) revealed local cr ystalline order consistent with a face-centered cubic structure.In the densest tapped configuration ν ∞ ＝0.675, a power-law singularity in g(r/d) near r/d＝1 was found, while at the lowest tapped solids fractionτ ∞ ＝0.583, the appearance of a double peak for 1.8<r/d<2.0resembled that seen in experiments and molecular dynamics simulations in the literature.The behavior of the mean coordination number, computed from the near-neighbor distributions as a function of the equilibrium bulk solids fraction was found to be in reasonably good agreement with experimental results in the literature.The ordering influence of the floor in the simulations was pronounced, but its occurrence is in accord with experimental observations in the literature.Preliminary calculations of the solids fraction behavior from the floor demonstrated the ordering effect of a smooth base.Further investigations of this effect are underway for deeper systems with a larger number of particles under various tapping strategies.The findings of this study with regard to the effect of tap intensity on densification support our previous conjecture 56) on the existence of a range of vibration parameters that favors an increase of solids fraction in continuously oscillated systems.Our Monte Carlo results suggest that a bulk solid will experience an increase in density for a range of tap intensities that promote a 'slow' rearrangement of particle positions, i.e., for tap intensities that do not inhibit the creation of local ordering in the microstructure.On the other hand, for larger accelerations for which the system is aggressively expanded at each tap, the bulk solids fraction does not grow as a function of the number of taps.Discrete element studies are currently underway to study these phenomena.

Fig. 3
Fig. 3 Standard deviation of the bulk solids fraction〈σ(n;γ)〉versus tap number N T at tap intensity γ＝1.25 for 100 realizations.This demonstrates the convergence properties of the methodology.

Fig. 4 Fig. 2
Fig.4 Dependence of the ensembled-averaged bulk solids fraction on the number of realizations M at tap intensity γ＝1.25 to exhibit the sensitivity of the system to the ensemble size.The vertical lines through the points are error bars.

Fig. 8
Fig.8 Comparison of the evolution of the ensemble-averaged bulk solids fraction at intensity γ＝1.25 with(×) and without(□) the application of random lateral displacements.The results with and without these displacements are statistically indistinguishable.

12Fig 13
Fig 13 Solids fraction〈ν〉as a function of normalized height y/d from the floor for the case ν ∞ ＝0.654.

Table 1 .
Fitting parameters for hyperbolic tangent law

Table 2 .
Fitting parameters for stretched exponential law

Table 3 .
Fitting parameters for inverse logarithmic law

Table 4 .
Fitting parameters for reciprocal linear law