Experiments and Simulations on Vibration Induced Densification of Bulk Solids

An important industrial issue in the production, transport and packaging of granular materials is the identification of those factors that affect bulk density, including how the material has been handled throughout its processing history. A common phenomenon that is not well-understood is the increase in bulk density that takes place when a container of granular materials is shaken, tapped or vibrated. There is a rather extensive literature on this topic, ranging from very fundamental investigations on the development of packing algorithms and microstructural analysis to more practical experimental and numerical studies. While it is not possible in this paper to provide a complete review of the extensive literature on this topic, a few representative papers are discussed in an attempt to provide a rough sketch. From a certain perspective, one may think of the phenomenon of density relaxation as having its foundations in the extensive literature on the packing of spheres. (See for example reference1-15). As early as 1611, Kepler explored the geometry of the snowflake while in 1665, Robert Hook investigated circle and sphere packings. In 1694, Gregory, a Scottish astronomer, suggested that thirteen rigid uniform spheres could be packed around a sphere of the same size a hypothesis that was only disproved by Leech some 262 years later. In 1727, Hales examined the packing of dry peas pressed into a container an experiment that was later known as the ‘peas of Buffon’ based on work in 1753 by Comte de Buffon. In 1887, Thompson considered the question of how to fill Euclidean space using truncated octahedrons, while Slichter was the first to attempt to find analytical expressions for the porosity in beds of uniform spheres. There has also been interest in the alternative question of what is the minimal solids fraction vm of a rigid assembly of uniform spheres such that each must touch at least four others and the contact points must not lie all in one plane. Along these lines, in 1932 Hilbert found a structure for which, vm=0.123, while a year later, Heesch and Laves created a stable arrangement of spheres such that vm=0.056. It should Abstract


INTRODUCTION
An important industrial issue in the production, transport and packaging of granular materials is the identification of those factors that affect bulk density, including how the material has been handled throughout its processing histor y. A common phenomenon that is not well-understood is the increase in bulk density that takes place when a container of granular materials is shaken, tapped or vibrated. There is a rather extensive literature on this topic, ranging from very fundamental investigations on the development of packing algorithms and microstructural analysis to more practical experimental and numerical studies. While it is not possible in this paper to provide a complete review of the extensive literature on this topic, a few representative papers are discussed in an attempt to provide a rough sketch. From a certain perspective, one may think of the phenomenon of density relaxation as having its foundations in the extensive literature on the packing of spheres. (See for example reference1-15). As early as 1611, Kepler explored the geometry of the snowflake while in 1665, Robert Hook investigated circle and sphere packings. In 1694, Gregory, a Scottish astronomer, suggested that thirteen rigid uniform spheres could be packed around a sphere of the same sizea hypothesis that was only disproved by Leech some 262 years later. In 1727, Hales examined the packing of dry peas pressed into a container -an experiment that was later known as the 'peas of Buffon' based on work in 1753 by Comte de Buffon. In 1887, Thompson 16) considered the question of how to fill Euclidean space using truncated octahedrons, while Slichter 17) was the first to attempt to find analytical expressions for the porosity in beds of uniform spheres. There has also been interest in the alternative question of what is the minimal solids fraction vm of a rigid assembly of uniform spheres such that each must touch at least four others and the contact points must not lie all in one plane. Along these lines, in 1932 Hilbert 18) found a structure for which, vm=0.123, while a year later, Heesch and Laves 19) created a stable arrangement of spheres such that vm=0.056. It should be pointed out, however, that these geometric structures cannot be created using standard experimental procedures. The distinction between a 'loose' (or poured) and 'dense' random packing of spheres was conceived by Oman and Watson 20) in 1944. This idea was quantified in experiments by Scott 21) , in which 3mm steel ball bearings poured into cylindrical containers were subjected to two minutes of shaking. Although the vibration parameters used were not repor ted, he found two values of the solids fraction, i.e., vloose=0.59 and vdense=0.63 corresponding to a random close and dense configuration. Improved experiments by Scott and Kilgour 22) yielded a more precise value of the dense solids fraction vdense=0.6366±0.0005. We note that deposition experiments in viscous liquids to achieve stable loose structures have reported solids fraction values as low as vloose=0.555±0.0005 for glass spheres and vloose=0.506±0.0005 for roughened acrylic spheres 23) . The notion that a consolidated state of optimal bulk density could be achieved through the use of high frequencies and relatively small displacement amplitudes was suggested by Stewart in 1951 -a claim that found support in the 1964 experiments of Evans and Millman 24) . Several years earlier, Macrae et al. 25) proposed that bulk density was related to impact velocity, with a critical value producing optimal results. In 1967, D'Appolonia 26) per formed vibration tests on dry sand within a cylindrical vessel. A mechanical shaker was employed to produce unidirectional harmonic motion of the cylinder with displacement amplitudes up to 0.254 mm and frequencies 10 ＜ f ＜ 60 Hz. By measuring the volume change of the sand over a range of vibration parameters, a plot of bulk density versus dimensionless acceleration Г≡aω 2 /g was generated, indicating that the greatest increase took place at Г ～ ～2. It was unclear from their data why increasing the acceleration beyond an optimal value resulted in decrease in bulk density. Dobr y and Whitman 27) , who extended D'Appolonia's experiments, reported in 1973 that the most rapid compaction occurred when 0.9 ＜Г＜1.1, while a maximum bulk density was achieved for 1.1 ＜Г＜1.3. At higher accelerations between 1.3 and 2.0, the density either stabilized or continued to increase. In addition to the use of continuous vibrations, there have been investigations on tapped systems, such as that by Takahashi and Suzuki 28) who studied the evolution of the volume of real powders. They described the phenomenology via a first-order rate law whose solution yielded a solids fraction v(n) as a function of the number of taps n that evolved to an apparent final density v ∞ in accordance to ν(n)=ν0+(ν∞-ν0)e -k/n . The effect of detached, vertical sinusoidal taps applied to a tall cylindrical vessel filled with 2mm mono-disperse, soda-lime glass spheres was experimentally studied by Knight et al. 29) using a noninvasive, capacitive technique to measure solids fraction. The evolution of the measured solids fraction was found to rely on the relative acceleration and on the number of taps. Experimental data was fit to a four-parameter phenomenological model of the form ν(n)=ν ∞-ν ∞-ν 0 1+Bln(1+n/τ) where v ∞ is the steady-state density (dependent on the acceleration history), τ is a relaxation time, and B is an undetermined constant that depends on Г. Linz 30) proposed an explanation of the experimental model from an analysis of the stroboscopic decay law that he derived from his physical interpretation of the compaction process. Knight et al.'s experiments were extended by Nowak et al. 31,32) to explore the frequency dependence and amplitude of the measured density fluctuations as a function of vibration intensity Г. They found that at certain intensities the system attained a well-defined average steady-state density with large fluctuations after extended tapping. The magnitude of these fluctuations depended not only on the depth at which measurements were made, but also on Г (eg., increasing Г produced larger fluctuations about the mean density). In addition to physical experiments, particle level simulations that offer insights into the micro-structural process taking place during the density relaxation process have been done. For example, Baker and Mehta 33) studied the qualitative effects of vibrations on a system of mono-disperse spheres using a hybrid simulation technique in which a real 'shake' was modeled by a controlled volume expansion of the assembly, a hard-sphere Monte Carlo at a low temperature to reduce the system potential energy, followed by a modified sequential random-close-packing algorithm to achieve stability. The authors carried out a careful analysis of the solids fraction and coordination number distribution as a function of their simulated shaking intensity, as well as a study of contact networks that corroborated the roles of various identified relaxation mechanisms. Rosato and Yacoub 34) carried out discrete element simulations to assess the effect low amplitude (a/d＜0.1), vertical oscillations applied to a vessel of frictional, inelastic spheres of diameter d. They reported fairly good fits of the data with phenomenological predictions of Knight et al. 29) and the exponential decay model of Takahashi and Suzuki 28) .
More recently, An et al. 35) presented discrete element results that supported the earlier findings of Zhang et al. 36,37) with regard to the effect of frequency and amplitude on solids fraction. They identified two respective densification mechanisms corresponding to low and high relative accelerations.
The physical system of interest in this paper is a model granular material comprised of acrylic spheres housed in a vertically oscillated cylindrical container. We propose a basis for trends obser ved in our experiments (Section 2) on the bulk solids fraction as a function of the applied vibration amplitude and frequency. Our explanation hinges on a detailed series of discrete element simulations (Sections 3 and 4) that qualitatively reproduce the experimental behavior. While the parameter space examined in the simulations is rather extensive (e.g., see reference 36)), in this paper we report on a subset of our studies that is relevant to the experiments. A phase diagram portraying 'improvement' in bulk density reveals distinct regions in the frequency-amplitude space which correlate with the experiments.

Experiments and Results
The containment vessel is an acr ylic cylinder of diameter D formed from stacked rings that are mounted onto a B&K shaker. Fig. 1 shows a schematic of the experimental system. An accelerometer attached to the piston provides feedback control through which precise adjustments to the frequency and amplitude could be made. The first part of the experiment involved measurements of the initial (before vibrations are applied to the piston) bulk solids fraction. Particles (acr ylic spheres, d = 3.175 mm) were slowly poured into the cylinder and then the top layer was removed by sliding the top ring across. This was done to ensure that the cylinder was filled with a level surface to a height of 30d. The poured bulk solids fraction ν 0 could then be computed from the volume of the cylinder and weight of its contents. For the aspect ratio D/d=20 used in the experiments, we obtained an average value ν 0 = ～ 0.604 that is typically associated with a loose random packing 6) .
The vibration experiments were carried out by filling the cylinder with mono-disperse acrylic spheres (d = 3.175 mm) to an undisturbed bed depth H ～ ～ 95.3 mm using the method previously described. The filled cylinder was rigidly mounted onto the shaker head, which was sinusoidally oscillated over a range of amplitudes 0.04 ＜ a/d ＜ 0.24 and frequencies ω between 25Hz-100Hz, (corresponding to relative accelerations Г≡aω 2 /g between 0.94 and 11.54). For each selected frequency and amplitude, the shaker was run for ten minutes, with each experiment repeated several times to confirm the results. In order to reduce the buildup of static charge, the inside tube wall and particles were treated with a household antistatic agent. In measuring the bulk solids fraction ν1 after the vibrations were stopped, a level particle surface was formed using the same procedure as described in the pouring experiments. Figs. 2-6 show the improvement in solids fraction, ( ν1-ν0 ν1 ×100) versus Г for increasing values of a/d. Our experiments reveal four trends in the behavior of ν versus Г that depend on the level of the displacement amplitude. Fig. 2 shows that for a/d = 0.04, the solids fraction increases with Г, but the improvement is only about 3.3% at Г = 5.1. Observations of the vibrated bed indicated that there was little or no bulk motion with the exception of some activity of the particles at the top surface for Г＞ 2. We conjecture that at these excitation levels the bulk density increases through slow, cooperative particle rearrangements throughout the assembly, analogous to the 'push filling' mechanism suggested by An et al. 35) . When a/d is between 0.06 and 0.10, results in Fig. 3a show a maximum in bulk solids fraction (ν m ～ ～ 0.636, cor- responding to an improvement ＞ 5%) when Г is between 5 and 7, followed by a slight expansion of the bed (or a reduction in ν) with a further increase of Г. Near the peaks (at Г=Гc), particles near the cylinder walls were seen to remain in contact with each other and to travel downwards in a collective manner, a phenomenon indicative of convective motion that has been reported in other studies [38][39][40][41][42][43][44][45][46] . The movement of these wall particles became more pronounced as Г was increased to Гc. Further observations made near the walls and on the surface at Гc revealed closedpacked arrangements of the particles as illustrated in Fig. 3b. We believe that the influence of the cylinder walls in the presence of the slow convection in this region at Гc is responsible for the creation of these structures, which in turn contribute to the overall increase in bulk solids fraction. It is worthwhile noting that Nowak et al. 32) also pointed out that the walls of the cylindrical container used in their experiments (aspect ratio D/d ～ 9.4) may have influenced the compaction process. In fact, they attained a substantially larger mean solids fraction (ν= 0.656) than what is known as random close packing of spheres (v = 0.6366 6,22) ). Distinct oscillations (Fig. 4) in the solids fraction versus Г were found when a/d = 0.16, where the improvement in the bulk density was generally less than 3.6%. At this amplitude level, convection could still be observed as particles adjacent to the walls moved downward. At the highest amplitude level used (a/d = 0.24), the solids fraction versus Г remains almost constant and little improvement (approximately 0.5%) was attained (Fig. 5). At these conditions, a marked expansion of the bed depth occurred (analogous to other experiments reported in the literature 47) ) and the system resembled a dense, energetic "gas" of  particles (it required relatively little effort to insert a rigid bar into the vibrating mass). We were unable to see any particle convection near the cylinder walls in this case.

Discrete Element Simulations
The discrete element simulations reported in this paper employ the soft "partially latching spring model" developed by Walton et al. [48][49][50] for elastic-plastic collisions. This model features mechanisms of energy dissipation for both the normal (i.e., along the line of centers of the pair of colliding particles) and tangential directions. Energy loss in the normal direction is accomplished using linear springs that are activated as particles overlap and then move apart. Because the loading spring constant K1 is smaller than the unloading (or restituting) spring of constant K2, the normal relative separation velocity is smaller than the relative approach velocity, and this produces a constant effective restitution coefficient e = K1/K2 . For the flows simulated in this study, collision velocities were not of sufficient magnitude to warrant the use of a velocity-dependent restitution that is appropriate for particle velocities of the order of 1 m ／ s and greater.
In the tangential direction the Walton-Braun model approximates Mindlin's and Deresiewicz's theory 51) for elastic spheres subjected to tangential loading. Dissipation is achieved through a tangential stiffness that decreases with tangential displacement until it is zero, at which point full sliding occurs at the friction limit μ. Although contacts that experience rotation coupled with tangential sliding are not a feature of this model, particles can rotate due to the transmis-sion of tangential impulse. The time step Δt through which the particle equations of motion are integrated is approximated from the normal force model by dividing twice the time spent in unloading period during a particle collision into n steps. It can be shown that Δ t = πe n m /2K1 where m is the particle mass and e is the restitution coefficient. The value selected for K1 ensured that overlaps were less that approximately 0.01d in accordance with the behavior of real colliding particles, so that Δt ～10 -6 s. The mass density of the simulated particles (ρ=1200 kg/m 3 ) corresponded to acr ylic plastic to match the experimental material. The equations of motion of the N particles are integrated using a leap-frog method with a backward Euler approximation at t = 0. For the translational motion (rotation equations are analogous), the discretization is given by where F is the net force on the i th particle.
In the results of our study described next, the computational cell was a rectangular box having a solid side walls and a floor whose motion was governed by s(t)=asin(2π ft).

Simulation Results and Discussion
A poured assembly of spheres was obtained by  initially positioning particles randomly within the computational cell, after which they were allowed to fall under the action of gravity for 1.5 seconds. This generated a stable configuration with a bulk solids fraction that did not change appreciably for runs of longer duration. We remark that the results of our pouring simulations 36) produced bulk solid fractions that were consistent with experimental measurements in the literature. The poured assembly was then used as the star ting point for the vibration simulations. The assembly was energized by applying vertical, sinusoidal oscillations to the floor of the computational cell for three seconds, followed by a relaxation phase during which particles fell under gravity to a stable state. Although not reported in this paper, we found that the three seconds of vibration was sufficient for the system to reach steadystate conditions from time-averaged depth profiles of 'granular temperature' and solids fraction.
In what follows, the behavior of the bulk solids fraction as a function of amplitude and frequency is presented. In all cases, spheres were assigned a normal restitution coefficient e = 0.9, and friction coefficient μ p ＝0.1. As will be seen, the results were in reasonable qualitative agreement with the trends observed in our previously described experiments. In Section 4.2, a much larger system whose aspect ratio more closely matches the experiments is considered.
Here, the simulated results are presented in the form of a chart that shows the improvement in density as a function of frequency and amplitude.

Vibration Frequency and Amplitude
Four amplitude ratios were considered (a/d = 0.02, 0.08, 0.24, 0.48) over frequencies ranging from 5Hz to 90Hz. Spheres were assigned a normal restitution coefficient e = 0.9, and friction coefficient μ p ＝0.1. These values were obtained from experiments of Louge 52) for acrylic spheres (i.e., μ p ＝0.096 ± 0.006, e=0.934 ± 0.009). Although we reduced the value of e somewhat to improve computational efficiency, this did not produce any significant changes in the qualitative trends reported in this paper. The selection of a shallow configuration (i.e., poured fill height of approximately 7d andν0 ～ ～ 0.577) and an aspect ratio L/d = 9.4 minimized the system size (N = 600), so that greater computational efficiency was obtained. Despite the small system size, the qualitative trends obtained were consistent with those found in the physical experiments. Although not reported here, we also carried out studies in deeper beds and found similar behavior of the system's 'dynamic' state 36) , which in turn impacts the densification upon relaxation.
The general trend in Fig. 6 when a/d = 0.02 is very much the same as what took place in the physical experiments (see Fig. 2). Over the frequency range tested at this amplitude, there is a continual improvement in bulk solids fraction. When a/d = 0.08, the data shown in Fig. 7 suggests a peak around 6% near f = ～ 50Hz; for greater frequencies, a reduction in solids fraction takes place. When the amplitude a/d is 0.24 (Fig. 8), the peak value occurs at approximately 40 Hz, and the cur ve decays thereafter until, near 80 Hz, no improvement in bulk density is possible at higher frequencies. The occurrence of the peak and decay afterwards is consistent with the experimental obser vations of Appolonia et al. 26) . Finally, at a/d= 0.48, the improvement in solids fraction is minimal (Fig. 9), and after a frequency of approximately 35 Hz, the system does not experience any densification upon relaxation. This trend is analogous to the experiments reported in Fig. 5.  A comparison between the experiments (Fig. 2-5) and the simulation (Fig. 6-9) shows reasonable qualitative agreement. Furthermore, the maximum improvement in solids franction of approximately 6% agrees with the experimental measurements.
Although there are quantitative differences between the simulated and experimental results (possibly attributed to boundary conditions and aspect ratio), the simulation does generate all of the important critical phenomena observed in the experiments.

Densification Phase Chart
We studied a system of N = 8000 spheres in a computational cell having an aspect ratio L/d ～ ～ 25. In so doing, we were able to further validate the simulations in a system whose cross-sectional dimension was comparable to our experiments. The procedure previously described was followed to generate the initial poured assembly (ν0 ～ ＝0.604), which filled the cell to a depth of approximately 11d. Due to limited computational resources, we were unable to carr y out studies for deeper beds at this aspect ratio, which would have permitted real quantitative comparisons with our experiments.
We ran a test to identify a specific value of the amplitude and frequency at which the relaxed bulk density was largest; this occurred at a/d = 0.16 and f = 40 Hz. With these parameters, the system was vibrated until the computed solids fraction curve flattened out (Fig. 10). An extrapolation of the date as 1/t → 0 yielded a solids fraction ν ＝0.658, in close agreement with the experimental results of Nowak et al. 53) . Arrangements of the particles adjacent to the side walls at t = 13s revealed the formation of hexagonally packed structures analogous to what was observed in  our experiments A series of case studies was carried out spanning vibration amplitudes and frequencies 0.02 ＜ a/d ＜ 0.48 and 10Hz ＜ f ＜ 90Hz. Our results are summarized in Fig. 11 in the form of a chart showing the improvement in bulk solids fraction as a function of amplitude and frequency. Four distinct regions are visible corresponding to various levels of improvement as indicated by the gray scale. A region of opti-mal improvement (in the amplitude-frequency space) is clearly visible, a trend that is in agreement with our experiments. Although not presented in this paper, we also examined profiles of the system's granular temperature, solids fraction and ratio of lateral to vertical kinetic energy, which revealed common characteristics for each region 37) . At the poorest levels of improvement, the vibrated system was found to be in a relatively high energetic state which settled to a  static configuration with little or no increase in bulk density. The full details of our investigation correlating the vibrated system with its relaxed density will be the subject of a future publication.

Conclusions
The densification behavior that takes place when a ver tically oscillated system of uniform spheres is allowed to relax under gravity was investigated.
The work finds its motivation in the industrial sector involved with the packaging of bulk solids. Experiments were carried out that exhibited clear trends regarding the influence of amplitude and frequency in achieving a state of optimal density. These trends were qualitatively reproduced in discrete element simulations that were extended further to a system having an aspect ratio comparable to the physical experiments. The results of the latter study were summarized in a 'phase' chart that revealed distinct frequency-amplitude regions which were characterized by various levels of the 'improvement' in the bulk density. Except for quantitative differences that depend on particle properties and mass overburden, we expect that similar patterns for the improvement would emerge for other materials.