Simulation of Inhomogeneous Columns of Beads under Vertical Vibration

Abstract

We present results from computational experiments of one-dimensional inhomogeneous granular column. With characteristic properties lying between those of liquids and solids, this class of material consists of assemblies of solid particles which interact via mechanical forces (contact and friction) and are maintained together by a gravitational field. The system is composed of a vertical column of spherical beads driven by a sinusoidally vibrating plate. Results for 10 bead experiments show the fluidization and condensation phenomena, both dependent on the amplitude and frequency of the driving plate. Some scaling properties arising from the dynamics of the beads are also shown.

1. Introduction

Granular materials are present in everyday life playing a critical role in several fields of applications such as pharmaceutical, cosmetic, and housing industries. Knowlegde of granular materials is essential to improve procedures involving extraction, transport, storage, and mixing processes. Each process can present desired and undesired phenomena occasioned by intrinsic features of the granular dynamic. Moreover, understanding these proprieties turns out to be important to minimize handling costs such as waste of energy. Granular flow systems can present clogging such that mixture procedures may produce segregation when handling the material in an inadequate manner [1]. Due to the large and practical applicability, this field was first restricted to engineering studies. However, as granular systems can also exhibit period doubling followed by chaotic behavior, they have being investigated by physicists and computational mathematicians.

Earlier studies [2–4] have investigated the collective behavior of a single column of beads composed by identical spheres, with the same radius and mass. These works rely on a simple model with minimum degree of details, namely, neglecting rotational movement and tangential components to describe specific features of the granular material. Results display different regimes similar to solids and fluids just by adjusting a single parameter. Other analyses quantify dilation of the column and energy transfer to the system, which are also important to model the behavior of the granular column.

Later works focused on describing the granular dynamic by using more realistic models. Some authors extended one-dimensional studies to two-dimensional columns and included rotational movement in the model [5, 6] and verified through a two-dimensional study and in the absence of gravity clustering formation using a velocity dependent coefficient of restitution [7]. Other improvements include applying more realistic interaction contact models taking into account different tangential forces such as the ones proposed by Cundall and Strack [8], Walton and Braun [9].

In this paper, we investigate the occurance of fluidization and condensation regimes in a column of beads having an inhomogeneous distribution of mass. In the fluidized regime, the motion of the individual particles looks erratic, and the beads are well separated, while in the condensed regime, the motion becomes more regular with the beads moving as a whole and locked on the excitation frequency of the drive source [2]. As is widely found in granular mixtures of particles with different densities of mass, a change in mass distribution significantly impacts the dynamics of the system upon the inclusion of an additional control parameter, namely, the beads mass ratio. Since clustering and condensation regimes of nearby particles consist of a collective motion with low or nearly zero relative velocities, changes in the velocity equations can modify the onset conditions of such processes.

The present paper is organized as follows. Section 2 explains both the model and the numerical procedures used to simulate a column of beads driven by a sinusoidally vibrating plate. Results are given in Section 3 which addresses discussions on fluidization and condensation states as well as scaling laws related to the center of mass of the column. Conclusions are presented in Section 4.

2. Numerical Procedure

The simulation system consists of a vertical column of spheres with radius 𝑟 and mass 𝑚_𝑖, 𝑖=1,2,…,𝑁 driven by a sinusoidally oscillating plate of infinitely large mass, as shown in Figure 1. Spheres are allowed to move only in the vertical direction and, therefore, they cannot interchange positions. The coefficient of restitution is assumed to be constant and, therefore, dissipation only occurs during collisions between particles with a coefficient of restitution 0<𝑒<1. Dissipation may also occur between the lowest sphere and the plate, according to coefficient of restitution 𝑒_𝑤.

Figure 1 

Model of a vertical column of spheres under gravity and driven by a vibrating plate.

Collisions between spheres are instantaneous and head-to-head since there is no angular momentum transfer between particles. The post collisions velocities follow the relations 𝑣₁=𝑚₂𝑢₂(1+𝑒)+𝑢₁(𝑚₁−𝑚₂𝑒)𝑚₁+𝑚₂,𝑣₂=𝑚₁𝑢₁(1+𝑒)+𝑢₂(𝑚₂−𝑚₁𝑒)𝑚₁+𝑚₂, (2.1) where 𝑢_𝑖 and 𝑣_𝑖 are the corresponding velocities of the 𝑖th particle before and after collision. A similar treatment quantifies the interaction of the lowest sphere with the base, which is driven by a sinusoidal signal 𝐴sin(𝜔𝑡), where 𝐴 is the amplitude and 𝜔=2𝜋𝑓 the angular frequency. The instant of collision 𝑡^∗ between the bottom sphere and the base has to be determined by solving a system of equations of motion doing 𝑧(𝑡^∗)=𝑆(𝑡^∗) as follows: 𝑧(𝑡)=𝑧₀+𝑣₀𝑡+𝑔(𝑡)²2,𝑆(𝑡)=𝐴sin(𝜔𝑡). (2.2)

This nonlinear model is often referred to as the bouncing ball model studied due to the richness of nonlinear and chaotic phenomena [10–12]. In order to avoid numerical problems in the computation, it is more convenient to use dimensionless variables. They are obtained according to the scheme:^0𝑥0𝑡^02𝑑𝑐≡𝜏=𝜔𝑡,𝑧^{0𝑥002𝑑𝑐}=𝑧𝑔/𝜔²,𝑣^{0𝑥002𝑑𝑐}=𝑣𝑔/𝜔,Γ=𝐴𝑤²𝑔, (2.3) such that (2.2) reduces to^0𝑥0𝑧^02𝑑𝑐(𝜏)=𝑧^{0𝑥002𝑑𝑐}₀+𝑣^{0𝑥002𝑑𝑐}₀𝜏𝜏−²2_{0𝑥002𝑑𝑐},𝑠(𝜏)=Γsin(𝜏), (2.4) where Γ is the nondimensionalized acceleration.

By using numerical methods [2], the next collision time 𝜏^∗ is calculated, and the postcollision velocity is updated according to following collision rule:^0𝑥002𝑑𝑣^𝑐_{00𝑥002𝑑𝑐}𝑣₁=101+𝑒_𝑤−𝑒_𝑤^{0𝑥002𝑑𝑐}_00𝑥𝑢_002𝑑𝑐𝑢₁, (2.5) where ^0𝑥0𝑢^02𝑑𝑐₀ and ^0𝑥0𝑣^02𝑑𝑐₀ are the velocities of the base, and 𝑒_𝑤 is the coefficient of restitution between the sphere and the base. The collision time is calculated through Newton's method, from which the crossings of a sinusoidal function with a parabolic curve are determined. The nondimensionalization procedure ensures finding only two roots, one negative and the other positive, where the later is the solution sought.

A similar procedure holds for collisions between the upper spheres, 2≤𝑖≤𝑁, which follow between events (a collision or a sudden change of momentum) the equation of motion^0𝑥0𝑧^02𝑑𝑐_𝑖(Δ𝜏)=𝑧^{0𝑥002𝑑𝑐}_𝑖,0+𝑣^{0𝑥002𝑑𝑐}_𝑖,01Δ𝜏−2(Δ𝜏)² (2.6) with Δ𝜏 denoting the time span between two consecutive collisions, the first of which at the time 𝜏₀ assigning the initial position 𝑧_𝑖,0 and initial velocity 𝑣_𝑖,0 of a generic sphere. The second event or collision between two spheres, is identified by ^0𝑥0𝑧^02𝑑𝑐_𝑖(Δ𝜏)=𝑧^{0𝑥002𝑑𝑐}_𝑖−1(Δ𝜏) which givesΔ𝜏=−^{0𝑥002𝑑𝑐}𝑧_𝑖,0−𝑧^{0𝑥002𝑑𝑐}_{𝑖−1,00𝑥002𝑑𝑐}𝑣_𝑖,0−𝑣^{0𝑥002𝑑𝑐}_𝑖−1,0, (2.7) where ^0𝑥0𝑧^02𝑑𝑐 and ^0𝑥0𝑣^02𝑑𝑐 are the position and the velocity of 𝑖th particle, respectively. A negative solution Δ𝜏<0 means a collision in the past (at a time prior to 𝜏₀) or, equivalently, that the two spheres are moving away from each other, while a solution Δ𝜏=0 implies that the spheres remain in contact. Solutions of interest are those Δ𝜏>0 thus meaning the occurance of a forward collision. Once Δ𝜏 is determined, the position and velocity of all spheres are updated, and the system is released thereafter to evolve according (2.6) until a next event is identified [2]. All calculated values are stored, and the smallest positive 𝜏^∗ determines the next collision time, which, in turn, is numerically iterated to calculate the initial values for the forward collision. In between events, we assume that all particles describe trajectories following the classical laws of mechanics. Positions and velocities of all particles must be updated in each event, and the ordered list of time events is managed by adding the latest event times.

3. Results

First, setting the parameters 𝑁=2 spheres, Γ=1.0, 𝑒_𝑤=1, 0.6<𝑒<0.9 (𝑒 is the restitution coefficient for collision between spheres), and 𝑚_𝑏/𝑚_𝑡=2 (indices 𝑏 and 𝑡 correspond, resp., to bottom and top) leads to a steady-state motions in which both of the balls execute 1-periodic motions which are out of phase by half a period. This is a nice example of a ball bouncing over another bouncing ball, which impacts on the harmonically driven table. In performing the simulation, the two balls are dropped from arbitrary heights and after the transient is finished, the balls synchronize mutually by executing steady-state motions that consist of equal jumps shifted by a half oscillation period.

This periodic mode is stable and presents a large basin of attraction for dissipation values, as we found periodic convergence within the interval of 0.6<𝑒<0.9. For matter of simplicity, we keep 𝑒_𝑤 unit to quantify the overall dissipation by using a single parameter. Obviously, the amount of dissipated energy decreases when the restitution coefficient increases within the range 0.6<𝑒<0.9. This 2-ball system, in which the balls stay apart during their flights, bears a similarity to the segregation process in a binary granular mixture, wherein two granular materials differing by a physical property (be it, e.g., density, shape, or size) exhibit an inexorable propensity for segregation. In fact, this result is not achievable when we set a unit mass ratio, 𝑚_𝑏/𝑚_𝑡=1. We note that convection streams play a key role in the segregation of a two-dimensional granular system [13, 14], while in our one-dimensional model (Figure 2) segregation is accounted for timing and synchronization effects.

Figure 2 

At Γ=1.0, 𝑒=0.9, 𝑒𝑤=0.95, 𝑚𝑏/𝑚𝑡=2, we obtain an extended solution for the 1-mode periodic motion of the bouncing ball model. The lower ball hits the base and the upper ball at the same phase of oscillation.

Using 𝑁=10 spheres of the identical masses, Γ=1.5, 𝑒=0.92, and 𝑒_𝑤=1, we obtain the condensation regime with collective motion of particles for which the center of mass also describes a periodic motion (Figure 3). All the beads remain clustered together and remain collectively in phase with the driving plate. In a fashion quite identical to that of an inelastic ball, the beads move in the form of a single block (but not necessarily in contact), and the corresponding frequency spectrum of the center-of-mass motion displays subharmonics of the driving source, as shown by the inset in Figure 3. This regime no longer appears if we introduce an inhomogeneous mass distribution; motion of the beads is erratic, much like particles in a gas, with the center of mass describing a chaotic trajectory. Results for 𝑚_𝑖=2, 𝑖=1,…,5 and 𝑚_𝑗=1, 𝑗=6,…,10 are given in Figure 4. Decreasing Γ makes the particles come to a permanent contact with the plate. It has been verified also that a double mass ratio does not lead to the condensation regime.

Figure 3 

At Γ=1.5, 𝑚𝑡=𝑚𝑏, 𝑒=0.92, and 𝑒𝑤=1.0, the homogeneous string of ten balls is locked in the condensation regime. The frequency spectrum of the center-of-mass motion shows subharmonic responses from the driving source.

Figure 4 

At Γ=1.5, 𝑒=0.92, and 𝑒𝑤=1.0 (the same as in Figure 3), but placing 5 heavier spheres (mass 𝑚𝑏) on the bottom and 5 spheres (mass 𝑚𝑡) on the top of column with mass ratio 𝑚𝑏/𝑚𝑡=2, the spheres follow chaotic trajectories, as indicated in frequency spectrum of the center-of-mass motion.

Increasing the mass ratio by setting 𝑚_𝑖=10, 𝑖=1,2,…,5, and 𝑚_𝑗=0.1, 𝑗=6,7,…,10 produces a condensation regime coexisting with a fluidized regime as shown in Figure 5. The mass ratio is sufficiently large so that we can consider the column as two isolated systems. The center of mass, however, is essentially determined by the heavier spheres, which together with the center of the mass follows a periodic motion. Yet, we note that the center of mass fails to describe the behavior of the system accurately because the frequency spectrum does not show up the coexistence of the two regimes.

Figure 5 

At Γ=1.5, 𝑒=0.92, and 𝑒𝑤=1.0 (the same as in Figure 4) and placing 5 heavier spheres (mass 𝑚𝑏) on the bottom and 5 spheres (mass 𝑚𝑡) on the top of the column with mass ratio 𝑚𝑏/𝑚𝑡=100 makes the condensation and fluidization regimes to coexist (an intermediate situation between Figures 3 and 4). The inset shows the frequency spectrum of the center-of-mass motion.

From the previously assembled column in Figure 5 and by decreasing the coefficient of restitution to increase energy dissipation rate, we can drive the upper layer into the condensation regime. We set a coefficient of restitution for each layer making 𝑒₁=0.92 for the heavier particles and 𝑒₂=0.82 for the lighter ones. This strategy is useful when adjusting energy in each layer without interfering with another.

With this column so assembled, Figure 6 shows two groups of particles in the condensed regime. This result extends the result in Figure 3 for two layers moving in the form of two blocks; each one with its own characteristic weight.

Figure 6 

At Γ=1.5, 𝑒=0.92, 𝑒𝑤=1.0, and mass ratio 𝑚𝑏/𝑚𝑡=100 (the same as in Figure 5) the two groups of particles enter in the condensed regime upon adjusting the restitution coefficients to 𝑒𝑏=0.92 (heavier spheres) and 𝑒𝑡=0.82 (lighter spheres).

Next, using the mass ratio 𝑚_𝑏/𝑚_𝑡=100 as in Figure 5, but with 𝑒=0.67, and 𝑒_𝑤=1, we increase the acceleration parameter to a sufficiently large value (Γ=10) to fluidize the lower portion of the column (with five heavy spheres on the lower part). Results are shown in Figure 7, where the lower half part turns out fluidized, while the upper part (five light spheres) is in the condensed state because fewer energy (after being attenuated and absorbed by the bottom portion) rests available to fluidize the upper part. The upper layer (the lighter one), however, is more likely to be affected by smaller values of the restitution coefficient. For instance, decreasing this parameter to 𝑒=0.5, the light spheres become tightly clustered, like a single inelastic object, yet the lower spheres still remain in the fluidized state.

Figure 7 

At Γ=10, 𝑒=0.67, 𝑒𝑤=1.0, and mass ratio 𝑚𝑏/𝑚𝑡=100, the lower portion of the column becomes fluidized, while the upper half part remains in the condensed regime.

In order to quantify the degree of dilation, we use a convenient parameter defined by [2]𝜆=⟨𝑧_𝑛−𝑧₁⟩𝐴Γ, (3.1) which averages over time the column's height.

Then, we assembled two different columns with a specific global mass distributed among 10 spheres, which are arrayed into distinct arrangements. With a total mass of 20 m.u., the first arrangement includes five light spheres (1 m.u.) placed on the lower portion of the column and five heavy spheres (3 m.u.) on the top. Simulations with light spheres on the upper portion (with heavy ones on the bottom) as well as arrangements with light and heavy spheres alternated were performed. Also compared with homogenous columns (10 identical spheres, 2 m.u. each) in Figure 8, results indicate different degrees of dilation for each arrangement. We see that columns with light particles on the top show larger dilation, as linear momentum transfer from heavy particles to light ones allows the upper particles to jump higher. Arrangements with alternated spheres present the same dilation as that of the uniform distribution. Differently from the homogeneous column (green dots), the array with alternated spheres does not exhibit sawtooth oscillations at low values of the acceleration parameter (Γ). Configurations with heavy spheres on top have the smallest dilations since the collision energy is attenuated as it propagates up, thus making it impossible to fluidize the upper portion of the column. In all cases, and consistent with previous results [2], the column dilation remains constant when increasing the shaking parameter Γ>8.

Figure 8 

Plotted against the acceleration parameter Γ, time-averaged dilations of a homogeneous, and three inhomogeneous columns, 20 m.u. heavy, assembled with two sets of five spheres of mass ratio 1:3.

Results for the heavier column of 30 m.u., assembled with 5 m.u. and 1 m.u. spheres, are given in Figure 9, showing again that for homogeneous and alternating-sphere columns, the degrees of dilation are nearly the same as those for the 20 m.u. column (Figure 8). We see in Figure 9 that the heavier column (30 m.u.) with light spheres on the top exhibits a larger dilation than that obtained with the 20 m.u. column (Figure 8), since as expected, the heavier spheres on the bottom of the 30 m.u. column absorb energy from the table which, in turn, is transferred to the lighter spheres hovering on the top. For all cases, the dilation values converge to a constant value at large acceleration (Γ>10), meaning that the length of the dilated (or energized) column is proportional to 𝐴²𝜔², which is the time-averaged square velocity of the driving table.

Figure 9 

Plotted against the acceleration parameter Γ, time-averaged dilations of a homogeneous, and three inhomogeneous columns, 30 m.u. heavy, assembled with two sets of five spheres of mass ratio 1:5.

The amount of energy absorbed by the column is quantified in Figure 10, which compares the up shift of the center of mass of five columns assembled with ten spheres. The mass distribution of inhomogeneous columns is labeled by the mass ratio 𝑚_𝑡∶𝑚_𝑏 indicating that five spheres of relative mass 𝑚_𝑏 are placed on the lower half part of the column, and the remaining five spheres (of mass 𝑚_𝑡) are stacked on the upper part of the column. Simulations have been performed at the 20 Hz frequency, so that the dimensionalized center-of-mass up shift is expressed as Δℎ(cm)=6.20×10ℎ^{−20𝑥002𝑑𝑐}. Lower elevations are associated with the 2:1 column (with heavier spheres on the top and when the mass ratio 𝑚_𝑡/𝑚_𝑏 is largest), while the converse holds for the 1:10 columns. Moreover, a distinctive feature of column with 𝑚_𝑡/𝑚_𝑏<1 relates to the existence of an inflection point on the corresponding curves, which turn out concave down toward reaching a stagnation zone at about the critical value of the reduced variable 𝑋. The physical process behind the coalescence of the curves is simply the occurrence of condensation at values 𝑋_𝑐≈3.0, meaning that the whole column never becomes completely fluidized [2]. We also note that at around 𝑋_𝑐≈3.0, the center-of-mass up shift is close to zero, thus reflecting that the motion of the spheres is collective. For still larger values of 𝑋, the dissipation becomes very high, which makes the column act like a compact block with dynamic properties of a single inelastic object. On the other hand, at lower 𝑋 values, the columns become fully fluidized, with each curve reaching a saturation value.

Figure 10 

Displayed using the reduced variable 𝑋=𝑁(1−𝑒), the center of mass elevation ℎcm−ℎcm,0 (from the center of mass at rest) of 10-sphere columns of the same height assembled with spheres of different mass ratios.

4. Conclusions

The present paper has given, to our knowledge, the first exploratory study of the dynamics of inhomogeneous granular columns under vertical vibration. Relying on numerical methods and computer simulations, a number of features have been revealed, namely, the occurrence of fluidized and condensed states (resp., characterized by the acceleration parameter Γ and by the square velocity 𝐴²𝜔² of the of the driving source) and the convergence of the dilation curves to a region defined by a critical value of the reduced variable 𝑁(1−𝑒)≈3.0. Settings of control parameters (normalized acceleration Γ, coefficient of elastic restitution 𝑒, mass ratio 𝑚_𝑏/𝑚_𝑡) have been determined so as to drive the system in specific states, either through tailoring the column into separate arrays to oscillate collectively in phase or through fluidizing the whole column. In what seems to be the first reported example, a two-ball bouncing system has been discussed in connection with its control parameters. Driven in the 1-periodic oscillation mode, the system is quite stable over the range 0.6<𝑒<0.9, thus deserving a thorough analytical study by extension of the methods used in its one-ball counterpart [12, 15], upon including five equations in the dynamical mapping, in addition to a further control parameter, namely, the mass ratio 𝑚_𝑏/𝑚_𝑡 that arises from the mass inhomogeneity of the column.

Acknowledgments

This work has been supported by the National Council for Scientific and Technological Development (CNPq), Brazil. M. V. Carneiro also thanks the Physics Group of Federal University of Lavras (UFLA), MG, Brazil for computing facilities.