# Simulation of Inhomogeneous Columns of  Beads under Vertical Vibration

Authors: Marcus V. Carneiro; Joaquim J. Barroso; Elbert E. N. Macau
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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.


---

## Body

## 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. Section2
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 radiusr and mass mi, i=1,2,…,N 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<e<1. Dissipation may also occur between the lowest sphere
and the plate, according to coefficient of restitution ew.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(2.1)v1=m2u2(1+e)+u1(m1−m2e)m1+m2,v2=m1u1(1+e)+u2(m2−m1e)m1+m2, where ui and vi are the
corresponding velocities of the ith 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 Asin(ωt), where A is the
amplitude and ω=2πf the angular
frequency. The instant of collision t* between the
bottom sphere and the base has to be determined by solving a system of
equations of motion doing z(t*)=S(t*) as follows:  (2.2)z(t)=z0+v0t+g(t)22,S(t)=Asin(ωt).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:(2.3)t˜≡τ=ωt,z˜=zg/ω2,v˜=vg/ω,Γ=Aw2g, such that (2.2) reduces to(2.4)z˜(τ)=z˜0+v˜0τ−τ22,s˜(τ)=Γsin(τ), 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:(2.5)(v˜0v˜1)=(101+ew−ew)(u˜0u˜1), where u˜0 and v˜0 are the velocities
of the base, and ew 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≤i≤N, which follow between events (a collision or a sudden
change of momentum) the equation of motion(2.6)z˜i(Δτ)=z˜i,0+v˜i,0Δτ−12(Δτ)2 with Δτ denoting the
time span between two consecutive collisions, the first of which at the time τ0 assigning the
initial position zi,0 and initial
velocity vi,0 of a generic
sphere. The second event or collision between two spheres, is identified by z˜i(Δτ)=z˜i−1(Δτ) which
gives(2.7)Δτ=−z˜i,0−z˜i−1,0v˜i,0−v˜i−1,0, where z˜ and v˜ are the
position and the velocity of ith particle,
respectively. A negative solution Δτ<0 means a
collision in the past (at a time prior to τ0) 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 parametersN=2 spheres, Γ=1.0, ew=1, 0.6<e<0.9 (e  is
the restitution coefficient for collision between spheres), and mb/mt=2 (indices b and t  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 of0.6<e<0.9. For matter of simplicity, we keep ew 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<e<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, mb/mt=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, e=0.9, ew=0.95, mb/mt=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.UsingN=10 spheres of the
identical masses, Γ=1.5, e=0.92, and ew=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 mi=2, i=1,…,5 and mj=1, j=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, mt=mb, e=0.92, and ew=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, e=0.92, and ew=1.0 (the same as in
Figure 3), but placing 5 heavier spheres (mass  mb)  on the bottom
and 5 spheres (mass mt) on the top of
column with mass ratio mb/mt=2,  the spheres follow
chaotic trajectories, as indicated in frequency spectrum of the center-of-mass
motion.Increasing the mass ratio by settingmi=10, i=1,2,…,5,  and mj=0.1, j=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, e=0.92, and ew=1.0 (the same as in
Figure 4) and placing 5 heavier spheres (mass mb) on the bottom
and 5 spheres (mass mt) on the top of
the column with mass ratio mb/mt=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 Figure5 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 e1=0.92 for the heavier
particles and e2=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, Figure6 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, e=0.92, ew=1.0, and mass ratio mb/mt=100  (the same as in
Figure 5) the two groups of particles enter in the condensed regime upon adjusting
the restitution coefficients to eb=0.92 (heavier spheres)
and et=0.82  (lighter
spheres).Next, using the mass ratiomb/mt=100 as in Figure 5,
but with e=0.67,  and ew=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 e=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, e=0.67, ew=1.0,  and mass ratio mb/mt=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)λ=〈zn−z1〉AΓ, 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 of20 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 Figure9, 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 A2ω2, 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 Figure10, 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 mt:mb  indicating that
five spheres of relative mass mb are placed on
the lower half part of the column, and the remaining five spheres (of mass mt) 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 Δh(cm)=6.20×10−2h˜. Lower elevations are associated with the 2:1 column (with
heavier spheres on the top and when the mass ratio mt/mb  is largest),
while the converse holds for the 1:10   columns.  Moreover,
a distinctive feature of column with mt/mb<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 X. The physical process behind the coalescence of the
curves is simply the occurrence of condensation at values Xc≈3.0, meaning that the whole column never becomes
completely fluidized [2]. We also note that at around Xc≈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 X, 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 X values, the
columns become fully   fluidized, with each curve reaching a saturation value.Figure 10
Displayed using the reduced variableX=N(1−e), the center of mass elevation hcm−hcm,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 A2ω2 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 N(1−e)≈3.0. Settings of control parameters (normalized
acceleration Γ, coefficient of elastic
restitution e, mass ratio mb/mt) 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<e<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 mb/mt that arises
from the mass inhomogeneity of the column.



---

*Source: 345947-2009-04-22.xml*