Papers in Physics, vol. 7, art. 070007 (2015)

Received: 9 December 2014, Accepted: 13 April 2015
Edited by: L. A. Pugnaloni
Reviewed by: F. Vivanco, Dpto. de F́ısica, Universidad de Santiago de Chile, Chile.
Licence: Creative Commons Attribution 3.0
DOI: http://dx.doi.org/10.4279/PIP.070007

www.papersinphysics.org

ISSN 1852-4249

Density distribution of particles upon jamming after an avalanche in a
2D silo

R. O. Uñac,1∗ J. L. Sales,2 M. V. Gargiulo,2 A. M. Vidales1†

We present a complete analysis of the density distribution of particles in a two dimensional
silo after discharge. Simulations through a pseudo-dynamic algorithm are performed for
filling and subsequent discharge of a plane silo. Particles are monosized hard disks de-
posited in the container and subjected to a tapping process for compaction. Then, a hole
of a given size is open at the bottom of the silo and the discharge is triggered. After a
clogging at the opening is produced, and equilibrium is restored, the final distribution of
the remaining particles at the silo is analyzed by dividing the space into cells with different
geometrical arrangements to visualize the way in which the density depression near the
opening is propagated throughout the system. The different behavior as a function of the
compaction degree is discussed.

I. Introduction

Numerous studies have investigated the flow of
granular materials (such as glass beads, aggregates
or minerals, among others) through hoppers of vari-
ous geometries [1–8]. Commonly, the silo is initially
filled with the granular material at a given height.
Then, the outlet of the silo is opened and the mass
flow rate is recorded as a function of time. These
experiments provide useful information on the re-
lation between the flow rate and different geomet-
rical and physical properties of the system, such as

∗E-mail: runiac@unsl.edu.ar
†E-mail: avidales@unsl.edu.ar

1 Departamento de F́ısica, Instituto de F́ısica Aplicada
(UNSL-CONICET), Universidad Nacional de San Luis,
Ejército de los Andes 950, D5700HHW San Luis, Ar-
gentina.

2 Departamento de Geof́ısica y Astronomı́a. Facultad de
Ciencias Exactas F́ısicas y Naturales. Universidad Na-
cional de San Juan, Mitre 396 (E), J5402CWH San Juan,
Argentina.

the size and shape of the particles and the outlet,
the presence of friction forces and density fluctua-
tions [9]. Empirically, the flow rate is determined
by the known Beverloo’s equation which depends,
among other variables, on the apparent density in
the immediate neighborhood of the outlet region
of the silo. Furthermore, in the derivation of the
Beverloo’s equation, it is assumed that the region
primarily affected by the discharge is near the out-
let and has an effective diameter of the order of
the width of it, sometimes referred to as Beverloo’s
diameter [2].

Simulations have provided details about the
granular flow that are not accessible to experiments
such as the influence of frictional parameters be-
tween particles, particle shape, stress chains and
the statistics of the arches formed during a jam-
ming process [3, 4, 10–12].

When the size of the flowing particles is in the
order of the width of the aperture, jamming can
take place, and the particles stop flowing unless
additional energy is provided to the system. Many
practical problems are caused by jamming at the

070007-1



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

outlet of hoppers used in production lines, where
it is necessary to maintain a constant flow of ma-
terial. Analogous problems occur in the storage of
raw materials in silos, especially after a certain pe-
riod of accumulation or under manipulation opera-
tions that may cause a change in the packing frac-
tion. Thus, when it is necessary to empty the silo,
the material does not flow, either due to the pres-
ence of unwanted moisture or because of the com-
paction of grains sealing the outlet [13–17]. These
problems are caused by particle arches that form
at the outlet. Authors in Ref. [9] report data sug-
gesting temporal oscillations in the packing fraction
near the outlet. Those oscillations had a frequency
around 2 Hz. This result is important because it
is directly related to the likelihood of jamming in
the silo [4, 5]. Other authors performed flow ex-
periments using metal disks in a two dimensional
hopper in order to study the statistical properties
of the arches forming at the outlet [5]. Authors in
Ref. [7] found experimentally a linear variation for
the number of particles forming an arch with the
outlet size in a 2-D silo.

There are many works concerning the study of
packing density in flowing granular materials out
of a hopper. In particular, those dealing with the
presence of density waves and density fluctuations
in the bulk of the silo [18, 19]. Others focus on the
density distribution during discharge and analyze
the change of the density between the stagnant zone
and the flowing zone [20].

In a recent work [12], authors have studied the
jamming occurring in the flow through small aper-
tures for a column of granular disks via a pseudo-
dynamic model. The effect that the preparation
of the granular assembly has on the size of the
avalanches was investigated. To this end, pack-
ing ensembles with different mean packing frac-
tions were created by tapping the system at dif-
ferent intensities. This work succeeded in demon-
strating that, for a given outlet size, different mean
avalanche sizes are obtained for deposits with the
same mean packing fraction that were prepared
with very different tap intensities. Nevertheless,
a complete characterization of the density of par-
ticles inside the column, both before and after the
discharging process, was left out.

For all above, a study of the variation of den-
sity near the outlet of a silo, before and after the
discharge in the presence of a jamming, is impor-

tant to relate it to the subsequent behavior of the
system during the restart of the flow.

II. Simulation procedure

The implementation and description of a simula-
tion algorithm using a pseudo-dynamic code has
been developed in numerous studies since the lim-
inal work of Manna and Khakhar [12, 21–23]. As-
suming inelastic massless hard disks that will be
deposited in a 2D die simulating a silo, the pseudo-
dynamics will consist in small falls and rolls of the
grains until they come to rest by contacting other
particles or the system boundaries. We use a con-
tainer formed by a flat base and two flat vertical
walls. No periodic boundary conditions are applied.

The deposition algorithm consists in choosing a
disk in the system and allowing a free fall of length
δ if the disk has no supporting contacts, or a roll of
arc-length δ over its supporting disk if the disk has
one single supporting contact [12, 21, 22, 24]. Disks
with two supporting contacts are considered stable
and left in their positions. If during a fall of length
δ a disk collides with another disk (or the base), the
falling disk is put just in contact and this contact
is defined as its first supporting contact. Analo-
gously, if in the course of a roll of length δ, a disk
collides with another disk (or a wall), the rolling
disk is put just in contact. If the first support-
ing contact and the second contact are such that
the disk is in a stable position, the second contact
is defined as the second supporting contact ; other-
wise, the lowest of the two contacting particles is
taken as the first supporting contact of the rolling
disk and the second supporting contact is left un-
defined. If, during a roll, a particle reaches a lower
position than the supporting particle over which it
is rolling, its first supporting contact is left unde-
fined (in this way, the particle will fall vertically in
the next step instead of rolling underneath the first
contact). A moving disk can change the stability
state of other disks supported by it; therefore, this
information is updated after each move. The depo-
sition is over once each particle in the system has
both supporting contacts defined or is in contact
with the base (particles at the base are supported
by a single contact). Then, the coordinates of the
centers of the disks and the corresponding labels of
the two supporting particles, wall or base are saved

070007-2



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

for analysis.
An important point in these simulations is the

effect that the parameter δ has in the results since
particles do not move simultaneously but one at a
time. One might expect that in the limit δ → 0,
we should recover a fairly ”realistic” dynamics for
a fully inelastic non-slipping disk dragged down-
wards at constant velocity. This should represent
particles deposited in a viscous medium or carried
by a conveyor belt. We chose δ = 0.0062d (with d
the particle diameter) since we have observed that
for smaller values of δ, results are indistinguishable
from those obtained here [12, 24].

It is worth saying that the pseudo-dynamic algo-
rithm used here allows the final configurations ob-
tained after each tapping to be completely static,
given that each disk is supported by other two disks
as required by the equilibrium conditions in the
model. In this way, one can follow the history of
formation of the packing and, thus, the occurrence
of the arches. This leads to a straightforward defi-
nition of the arches in the system, as we will see be-
low. On the other hand, this algorithm is faster in
the generation of a given ensemble of particles and,
the subsequent tapping process, than the DEM one.
The ones above are the most important advantages
of the algorithm that justify our choice.

Because arch formation has been identified as a
potential cause for segregation in non-convecting
systems [25–27], we have centered previous research
on detecting arches and analyzing their behavior
and distribution in piles and jammed silos [23, 28].
Indeed, identification of arches is a rather complex
task and the results presented in those works were
novel and original at that time. We recommend to
those readers interested in the characterization of
the arches formed before and after the discharge of
a 2D silo filled with disks to address our previous
work in Ref. [12].

Arches are sets of mutually stabilizing particles
in a static granular sample. In our pseudo-dynamic
simulations we first identify all mutually stable par-
ticles and then find the arches as chains of parti-
cles connected through these mutual stability con-
tacts. Two disks A and B are said mutually stable
if A is the left supporting particle of B and B is
the right supporting particle of A, or vice versa.
Since the pseudo-dynamics rest on defining which
disk is a support for another disk during the de-
position, this information is available in our sim-

ulations. Chains of mutually stable particles can,
thus, be found straightforwardly. These chains can
have, in principle, any size starting from two par-
ticles. Details on the properties of arches found in
pseudodynamic simulations can be found in previ-
ous works [7, 24, 28].

III. Filling and emptying the silo

As explained above, the aim of the present work is
to analyze the density patterns of particles inside
a silo after its discharge and as a function of the
compaction degree before the avalanche event.

We first need to prepare packings at reproducible
initial packing fractions. To achieve this, we have
chosen a well known technique to generate repro-
ducible ensembles of packings [12]. Thus, we use a
simulated tapping protocol (see below) to generate
sets of initial configurations that have well defined
mean packing fractions. The simulations are car-
ried out in a rectangular box of width 24.78d con-
taining 1500 equal-sized disks of diameter d. Ini-
tially, disks are placed at random in the simula-
tion box (with no overlaps) and deposited using
the pseudo-dynamic algorithm. Once all the grains
come to rest, the system is expanded in the vertical
direction and randomly shaken to simulate a verti-
cal tap. Then, a new deposition cycle begins. After
many taps of given amplitude, the system achieves
a steady state where all characterizing parameters
fluctuate around equilibrium values independently
of the previous history of the granular bed. The
existence of such ”equilibrium” states has been pre-
viously reported in experiments [29].

The tapping of the system is simulated by mul-
tiplying the vertical coordinate of each particle by
a factor A (with A > 1). Then, the particles are
subjected to several (about 20) Monte Carlo loops
where positions are changed by displacing parti-
cles a random length ∆r uniformly distributed in
the range 0 < ∆r < A − 1. New configurations
that correspond to overlaps are rejected. This dis-
ordering phase is crucial to avoid particles falling
back again into the same positions. Moreover, the
upper limit for ∆r (i.e., A− 1) is deliberately cho-
sen such that a larger tap promotes larger random
changes in the particle positions. For each value
of A studied, 103 taps are carried out for equili-
bration followed by 5 × 103 taps for production.

070007-3



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

500 deposited configurations are stored which are
obtained by saving every ten taps during the pro-
duction run after equilibration. These deposits will
be used later as initial conditions for the discharge
and flow through an opening [12].

It is worth mentioning that friction is not consid-
ered in our present model, neither between particles
nor between particles and the walls (convection is
not present [12]). The key feature to stabilize the
particles inside the silo is the sticking to the base
of the bottom particles. This ensures the achieving
of a final static system. As one would expect, the
presence of friction would affect the packing den-
sity. In this sense, we have already studied the ef-
fect of a kind of adhesion force like capillary forces
in the packing fraction of tapped columns of disks
in a previous work [28].

For each deposit generated, we trigger a dis-
charge by opening an aperture of width D relative
to the diameter d of the disk in the center of the
containing box base. Grains will flow out of the box
following the pseudo-dynamics until a blocking arch
forms or until the entire system is discharged (with
the exception of two piles resting on each side of the
aperture). During the dynamics, disks that reach
the bottom and whose centers lie on the interval
that defines the opening will fall vertically (even
if the surface of the disk touches the edge of the
aperture). This prevents the formation of arches
with end disks sustained by the vertical edge of the
orifice. After each discharge, we record the final
arrangement of the grains left in the box.

One single discharge attempt is carried out for
each initial deposit. This allows us to assure that
the initial preparation of the pack belongs to the
ensemble of deposits corresponding to the steady
state of the particular tap intensity chosen. To il-
lustrate the initial and final packing configurations
inside the silo, we present two snapshots in Fig. 1
showing the particles before, part (a), and after,
part (b), the discharge for the case A = 1.1 and
D = 2.5.

In a previous work, we have focused on the effect
that the preparation of the granular assembly has
on the size of the avalanches obtained. In the next
section, we present the analysis of the correspond-
ing packing densities before and after the discharge.

  

(b)(a)

  

Figure 1: Snapshots of the packed particles inside
the plane silo, (a) before the discharge and (b) after
it. The case corresponds to A = 1.1 and D = 2.5.
Thin blue lines indicate contacts between particles
and thick red lines indicate the arches.

IV. Density sampling

To measure the density patterns, we perform the
analysis by dividing the packing space into cells,
following different geometrical criteria, and for A =
1.1, 1.5, 2.0, and D = 2.5. The first measure-
ments are performed over all the particles inside the
silo. We choose two particular arrangements for the
sampling cells, i.e., circular ring sectors centered at
the base, the maximum radius for the ring being
five times the width of the base, and radial angu-
lar sectors starting at the center of the base, where
the opening of the sector is ∆θ, with the inclina-
tion θ measured from the horizontal axis. Figure
2 (a) and (b) illustrate the cases. On the other
hand, and to especially focus on the region near
the outlet, we define in addition three different cell
configurations. They are shown in Fig. 2 (c)-(e).
The details are depicted in the figure caption.

In all the sampling cell analysis, the criterion for
density calculation was the same. The particles
whose centers fall into a given sector are counted,
and that number is then divided by the sector area,

070007-4



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

Figure 2: Sketch of the geometry of the cells used
to measure the density of particles inside the silo.
(a) Circular ring sectors of thickness ρ, centered at
the middle point of the base. (b) Radial angular
sectors starting at the center of the base, spanning
all the packing. (c) Sectors as those in part (a) but
delimited by two symmetrical lines at 60◦ respect to
the horizontal and a maximum radius of half width
of the base. (d) Sectors as in part (b) but limited by
a semicircle with radius of half width of the base.
(e) Column of rectangular sampling sectors, with
height ρ, and width equal to the outlet size.

giving thus the particle density in the sector. We
measure the density corresponding to the initial
packing and for the packing array after the dis-
charge for each one of the independent configura-
tions and obtain the mean values of the density over
those configurations presenting clogging. This was
repeated for for all the geometries used.

V. Results and discussion

In Fig. 3 we present the results for the particle den-
sity as a function of the distance of the ring to the
center (see Fig. 2 (a)). The thickness of each ring is
equal to 1.61 a.u., i.e., a particle diameter. In part
(a) of the figure we plot the behavior of the density
for the initial packings, before the discharge of the
silo is performed. As expected, a constant behavior
is observed as the radius increases. As also reported
elsewhere [28], the average packing density is larger
for A = 1.1 than for A = 1.5 or 2.0. The oscila-
tions observed for A = 1.1 are due to the presence
of order in the packing structure [24]. This order is
virtually absent for higher tapping intensities. The
decrease of the density at small distances from the
center, especially for A = 1.1 and 2.0, is a purely

0 20 40 60 80 100
0.0

0.1

0.2

0.3

0.4

0.5

0.1

0.2

0.3

0.4

0.5

(b)

 

 

 

R

 

 

 

 
D

is
ks

 d
en

si
ty

(a)

Figure 3: Particle density vs. the distance of the
ring to the center, using the cell depicted in Fig.
2(a). Here ρ = d. (a) Initial packing behavior. (b)
Final stage after discharge. A = 1.1 (red squares),
1.5 (blue up triangles) and 2.0 (green circles).

geometric effect related to the small area covered
by the smaller rings and the particular disposition
of the disks at the base. On the other hand, the de-
crease observed at large distances for A = 1.1 is be-
cause the initial arrangement of particles presents
a free surface which is tilted.

After the discharge (Fig. 3 (b)), the disk den-
sity steeply drops inside the region close to the
outlet, putting in evidence the size of the concave
hole formed after the avalanche. Besides, the den-
sity slightly falls with height, for all amplitudes.
To highlight this effect, we draw three lines corre-
sponding to the mean densities at the initial state
(Fig. 3 (a)). This effect is due to the lower density
of packing near the walls and it will be explained
later. Here again, the oscillations for A = 1.1 are
associated to the order in the packing structure.
A similar analysis can be done taking a different
thickness for the rings, giving results qualitatively
equal to the ones shown in Fig. 3.

In Fig. 4 we present the results for particle den-
sity as a function of the inclination angle of the sec-
tor respect to the horizontal for three values of A.
Part (a) corresponds to the initial packing structure
and part (b) shows the state after discharge. The
opening angle in each sector is 1◦. The horizontal
axis represents the angle of the most inclined side
of the sector, as indicated in Fig. 2 (b). The peaks
for A = 1.1 in both plots are associated with the

070007-5



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

0 20 40 60 80 100 120 140 160 180
0.2

0.3

0.4

0.5

0.6

0.3

0.4

0.5

0.6

(b)

 

 

 

  (º)

D
is

ks
 d

en
si

ty

 

 

 

(a)

Figure 4: Particle density vs. the inclination an-
gle of the sector respect to the horizontal, for three
values of A. (a) Initial packing behavior. (b) Fi-
nal stage after discharge. A = 1.1 (red squares),
1.5 (blue up triangles) and 2.0 (green circles). The
three horizontal lines indicate the initial average
densities for the three amplitudes.

ordered structure, showing important increments
of the density for 60◦ and 120◦ and, less impor-
tant, for 30◦ and 150◦. When increasing the an-
gular sector width, ∆θ, peak amplitude decreases,
virtually disappearing, keeping the qualitative be-
havior shown in Fig. 4. The average density for
A = 1.1 is slightly higher, as expected.

After discharge, a dip at the central part of the
curves appears as an evident consequence of the
created hole, and showing the presence of a cen-
tral region around the hole in which the density
is reduced. This reduction is more pronounced for
A = 1.1 because the mean avalanche size for that
tapping intensity is greater [12] and the affected
region is approximately between 60◦ and 120◦ for
all amplitudes. At the bottom part of the figure,
we plotted three lines indicating the initial average
densities for the three amplitudes to better visual-
ize the density change.

Although useful to have an overview for density
behavior, information is lost when averaging over
sectors spanning all packing configuration. For that
reason, and to analyze in more detail the region
close to the outlet of the bidimensional silo, we im-
plement three new arrangements for the sectors to
perform the density analysis, as indicated in Fig.
2 (c)-(e). Fig. 2 (c) shows a scheme for sectors

similar to those in part (a) but delimited by two
symmetrical lines at 60◦ respect to the horizontal
and a maximum radius of 20 a.u. (half-width of the
container). Part (d) sketches the same sectors as
in part (b) but limited by a semicircle with radius
20 a.u.. Finally, part (e) in the same figure, shows
a column of rectangular sampling sectors, each one
with height ρ, and width equal to the outlet size.

In Fig. 5 we present the results for A = 1.1 for
the density averaged over sectors as in Fig. 2 (d).
The upper part shows the density for the initial
packing structure as a function of the inclination
angle of the sector with ∆θ = 1◦ (squares) and
∆θ = 10◦ (circles). For comparison, the middle
part of the figure also shows the initial density but
for the sectors without the boundary semicircle. By
circumscribing the density calculation to a semicir-
cle centered at the outlet and whose radius is the
half-width of the system, the ordered structure of
the bottom part of the initial packing is more ev-
ident (peak occurrence) inside sectors from 60◦ to
120◦.

Analyzing the final state (bottom part of Fig. 5),
a visible decrease in density around the outlet and
also a large disorder is observed, especially between
60◦ and 120◦. Outside this range, the structure
of the packing seems virtually unchanged, thus re-
vealing the non-avalanche area. As before, the lines
indicate the initial average density for the packing.

In Fig. 6 we show the results for A = 1.5 with
∆θ = 1◦ (squares) and ∆θ = 10◦ (circles) . It is
clear from part (a) of the figure that density oscil-
lations are much smaller than those for A = 1.1,
even for ∆θ = 1◦. This proves the lack of order in
the packing structure, except for those sectors close
to the base of the packing.

After discharge, the density does not change sub-
stantially throughout the bulk, except for a slight
decrease in the vicinity of the outlet orifice. This
is shown in Fig. 6 (b) and it is more evident for
∆θ = 10◦. A slight density modulation can also be
observed.

As a partial conclusion, we can say that the den-
sity of the packing after discharge retains the look
of the original grain disposition, which is consistent
with the disorder obtained for that tapping ampli-
tude, i.e., lowering of the packing fraction [12]. The
results and conclusions for A = 2.0 are quite similar
to those for A = 1.5.

In Fig. 7 the results for density are plotted av-

070007-6



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

0.3

0.4

0.5

0.6

0 20 40 60 80 100 120 140 160 180
0.2

0.3

0.4

0.5

0.6

0.3

0.4

0.5

0.6

 

 

 

(a)

(c)
 

 

 

 

 (º)

D
is

ks
 d

en
si

ty (b) 

  

Figure 5: Density average vs. inclination angle for
A = 1.1. Squares correspond to ∆θ = 1◦ and circles
for ∆θ = 10◦. Upper: initial packing structure.
Middle: initial density but for the sectors without
the boundary semicircle, for comparison. Bottom:
final state after the discharge for the system in the
upper part. The horizontal line indicates the initial
average density for ∆θ = 10◦.

eraging over sectors as in Fig. 2 (c), for A = 1.1,
1.5 and 2.0, and ρ = 2 a.u.. As in previous cases,
the initial packing shows the periodicity related to
the ordered structure for small ρ (not shown here)
and small A. As ρ increases, fluctuations become
smaller and a practically constant value for pack-
ing density can be observed as a function of the
distance to the bottom center of the silo.

After discharge (Fig. 7 (b)), a sudden decrease
in density is observed up to a height equivalent
to the estimated Beverloo’s diameter, i.e., in our
present case, 4 a.u. Two regions can be distin-
guished. One corresponding to a radius of 2 a.u.,
where the density is zero, and the other, where the
density increases rapidly, almost reaching the ini-
tial bulk value.

For A = 1.5, the fluctuations are only present
for small R, near the base. After discharge, the

0 20 40 60 80 100 120 140 160 180
0.2

0.3

0.4

0.5

0.6

0.3

0.4

0.5

0.6

(b)

 

 

 

 

 (º)

  

 

D
is

ks
 d

en
si

ty

(a)

Figure 6: Density vs. inclination angle as in Fig. 5
(a) and (c), but for A = 1.5, ∆θ = 1◦ (squares) and
∆θ = 10◦. The horizontal line indicates the initial
average density for ∆θ = 10◦.

behavior is similar to the case A = 1.1, but with
the presence of much less fluctuations.

For all amplitudes, the depression in density is
confined inside the cone subtended by the lines at
60◦. Compare Fig. 7 with Fig. 4 to see that the
sectors corresponding to the stagnant zone keep the
initial disk density. A similar result is obtained for
A = 2.0.

Finally, we performed the density analysis on the
column made by sectors with a given thickness ρ, as
indicated in Fig. 2 (e). Figure 8 shows the results
for A = 1.1, 1.5 and 2.0, and ρ = 2 a.u.. In part (a)
of the figure, which corresponds to the initial pack-
ing structure, the density of disks remains constant
with height for different A. Fluctuations are again
related to the ordered structure and decreases for
greater A. In the case A = 1.5, as known, the initial
configuration is more disordered than for A = 1.1
(less fluctuations even for small ρ). The mean den-
sity of the packing coincides with the corresponding
values in previous figures.

070007-7



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

0 2 4 6 8 10 12 14 16 18 20
0.0

0.2

0.4

0.6

0.2

0.4

0.6

(b)
  

 

 

 

 

 

R

d
is

ks
 d

en
si

ty

(a)

Figure 7: Results for the density of particles aver-
aged over sectors as in Fig. 2 (c), for A = 1.1, 1.5
and 2.0, and ρ = 2 a.u.. (a) Initial packing. (b) Af-
ter the discharge. A = 1.1 (red squares), 1.5 (blue
up triangles) and 2.0 (green circles). The three hor-
izontal lines indicate the initial average densities for
the three amplitudes.

After the avalanche (Fig. 8 (b)), a change in den-
sity is observed up to a height of about 25 to 30 a.u..
Up to the order of 40 a.u., fluctuations decrease
significantly, while for greater heights (far from the
outlet) only a slight decrease of those fluctuations is
observed. This is probably related to the formation
of arches in the structure that prevents the occur-
rence of internal avalanches at greater heights. Fig-
ure 1 and Fig. 9 illustrate this point. Those figures
show the snapshots before and after the discharge
for A = 1.1 and A = 1.5, respectively. There, the
arches formed by the particles are indicated with
thick segments. In parts (b), a loose packing with
the presence of arches is the resulting structure af-
ter the avalanche [12].

The density drop in Fig. 8 (b) coincides with that
in Fig. 7 (b) and its extent is related to the region
mainly affected by the discharge, which is of the

0 10 20 30 40 50 60 70 80 90 100
0.0

0.2

0.4

0.6

0.2

0.4

0.6

(b)

  

 

 

 

 

 

H

d
is

ks
 d

en
si

ty

(a)

Figure 8: Density results on the column vs. height
H belonging to sectors with thickness ρ = 2 a.u.
(1.24d). A = 1.1 (red squares), 1.5 (blue up trian-
gles) and 2.0 (green circles). The three horizontal
lines indicate the initial average densities for the
three amplitudes.

order of one Beverloo’s diameter (4 a.u.). Besides,
a decrease of about 5% to 10% in the average bulk
density is observed, in coincidence with Fig. 3 (b).

The analysis of the column for A = 1.05 presents
some interesting features comparing to the preced-
ing results. Figure 10 shows the results for that
amplitude. First, we observe that the density fluc-
tuations vs. height in the arrangements (both ini-
tial and final) are higher, reinforcing the fact that
fluctuations increase with decreasing A. This is re-
lated with the higher order structure.

On the other hand, up to a height of 25 to 30 a.u.
(which involves not only the empty vault area but
the rarefaction caused by the avalanche), the den-
sity does not fluctuate for A = 1.5 and 2.0, while
it does fluctuate for A = 1.05 and 1.1. This is be-
cause these latter arrangements are more compact
and ordered, allowing propagation of the decreased
density with a certain periodicity (see Fig. 1). As

070007-8



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

( b )( a )(a) (b)

Figure 9: Snapshots of the particles inside the silo,
(a) before the discharge and (b) after it. The case
corresponds to A = 1.5 and D = 2.5. Thin blue
lines indicate contacts between particles and thick
red lines indicate the arches.

before, the radius of the empty vault coincides with
the Beverloo assumption.

It is noteworthy that the decrease in the average
density in the bulk after the discharge is only no-
ticed when analyzing the data in a column or in the
circular sectors of Fig. 2 (a), not in the other cases.
This would indicate that the main contribution to
the avalanche is done by the particles just above
the hole.

VI. Conclusions

According to the results shown so far, it can be said
that there is a clear area of rarefaction in the pack-
ing column after the avalanche discharge of the silo.
This area is focused on the outlet opening. This
lower density region has been analyzed with dif-

0 10 20 30 40 50 60 70 80 90 100
0.0

0.2

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

(b)

  

 

 

 

 

 

H

d
is

ks
 d

en
si

ty

(a)

Figure 10: Density results as in Fig. 8 but A =
1.05. The horizontal line indicates the initial aver-
age density.

ferent geometries that provide a structural macro-
scopic view of it.

The presence of arches before and after the dis-
charge allows relating the packing structure with
the size of the density depression [12].

The higher the order in the initial structure, the
greater the extent of the rarefaction area after dis-
charge (Fig. 7 (b) and Fig. 8 (b)), i.e., a more
open (less ordered) initial structure induces a less
spreading of the lower density region.

Regarding the shape of the rarefaction region,
the decrease in density is more pronounced upward
and sideways to an angle of 60◦ (120◦) (which de-
fines the stagnation zone). This is evidenced by
comparing the density profiles for the columns (Fig.
8 (b)) and those for the circular sectors limited by
straight lines at 60◦ and 120◦ (Fig. 7 (b)) with
respect to the case shown in Fig. 5 (c).

After discharge in the case of Fig. 7 (b), a sudden
decrease in density is observed up to a height equiv-
alent to the estimated Beverloo’s diameter. Two
regions can be distinguished: one corresponding to

070007-9



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

a radius of 2 a.u., where the density is zero, and
the other where the density increases rapidly, al-
most reaching the initial bulk value.

The results for the case D = 2.75 are qualita-
tively the same as the one analyzed above, for that
reason they are not presented here.

It is important to remember that here we con-
sider a monosized distribution. Taking into account
that in real applications the size distribution of par-
ticles is usually not monodisperse, it is a future
challenge to see how our present results are modi-
fied when considering a given dispersion in the size
distribution of particles.

Acknowledgements - This work was supported
by CONICET (Argentina) through Grant PIP 353
and by the Secretary of Science and Technology
of Universidad Nacional de San Luis, Grant P-3-1-
0114.

[1] A W Jenike, Gravity Flow of Bulk Solids,
University of Utah Engineering Experiment
Station, Bull. 108 (1961), and Bulletin 123
(1964).

[2] W A Beverloo, H A Leniger, J van de Velde,
The flow of granular solids through orifices,
Chem. Eng. Sci. 15, 260 (1961).

[3] J Wu, J Binbo, J Chen, Y Yang, Multi-scale
study of particle flow in silos, Advanced Pow.
Tech. 20, 62 (2009).

[4] G H Ristow, Outflow rate and wall stress for
two-dimensional hoppers, Phys. A 235, 319
(1997).

[5] K To, P-Y Lai, Jamming pattern in a two-
dimensional hopper, Phys. Rev. E 66, 011308
(2002).

[6] S-C Yang, S-S Hsiau, The simulation and
experimental study of granular materials dis-
charged from a silo with the placement of in-
serts, Pow. Tech. 120, 244 (2001).

[7] A Garcimart́ın, I Zuriguel, L A Pugnaloni,
A Janda, Shape of jamming arches in two-
dimensional deposits of granular materials,
Phys. Rev. E 82, 031306 (2010).

[8] R O Uñac, O A Benegas, A M Vidales and I
Ippolito, Experimental study of discharge rate
fluctuations in a silo with different hopper ge-
ometries, Pow. Tech. 225, 214 (2012).

[9] R L Brown, J C Richards, Profile of flow
of granulates through apertures, Trans. Inst.
Chem. Eng. 38, 243 (1960).

[10] P W Cleary, M L Sawley, DEM modelling of
industrial granular flows: 3D case studies and
the effect of particle shape on hopper discharge,
Appl. Math. Model. 26, 89 (2002).

[11] A Anand, J S Curtis, C R Wassgren, B C Han-
cock, W R Ketterhagen, Predicting discharge
dynamics from a rectangular hopper using the
discrete element method (DEM), Chem. Eng.
Sci. 63, 5821 (2008).

[12] R O Uñac, A M Vidales and L. A. Pugnaloni,
The effect of packing fraction on the jamming
of granular flow through small apertures, J.
Stat. Mech. P4008 (2012).

[13] C Mankoc, A Janda, R Arévalo, J M Pastor,
I Zuriguel, A Garcimart́ın, D Maza, The flow
rate of granular materials through an orifice,
Gran. Matt. 9, 407 (2007).

[14] A P Huntington, N M Rooney, Chemical En-
gineering Tripos Part 2, Research Project Re-
port, University of Cambridge (1971).

[15] F C Franklin, L N Johanson, Flow of granular
material through a circular orifice, Chem. Eng.
Sci. 4, 119 (1955).

[16] S Humby, U Tüzün, A B Yu, Prediction of
hopper discharge rates of binary granular mix-
tures, Chem. Eng. Sci. 53, 483 (1998).

[17] A Mehta, G C Barker, J M Luck, Cooperativity
in sandpiles: statistics of bridge geometries, J.
Stat. Mech. P10014 (2004).

[18] G H Ristow, H J Herrmann, Density patterns
in two-dimensional hoppers, Phys. Rev. E 50,
R5 (1994).

[19] A Medina, J A Córdova, E Luna, C Treviño,
Velocity field measurements in granular gravity
flow in a near 2D silo, Phys. Lett. A 250 111
(1998).

070007-10



Papers in Physics, vol. 7, art. 070007 (2015) / R. O. Uñac et al.

[20] L Babout, K. Grudzien, E Maire, P J Withers,
Influence of wall roughness and packing den-
sity on stagnant zone formation during funnel
flow discharge from a silo: An X-ray imaging
study, Chem. Eng. Sci. 97, 210 (2013).

[21] SS Manna, D V Khakhar, Internal avalanches
in a granular medium, Phys. Rev. E 58, R6935
(1998).

[22] Manna S S, Self-Organization in a Granular
Medium by Internal Avalanches, Phase Tran-
sit. 75, 529 (2002).

[23] R O Uñac, J G Benito, A M Vidales, L A Pug-
naloni, Arching during the segregation of two-
dimensional tapped granular systems: Mix-
tures versus intruders, Eur. Phys. J.E 37, 117
(2014).

[24] L A Pugnaloni, M G Valluzzi and L G Valluzzi,
Arching in tapped deposits of hard disks, Phys.
Rev. E 73, 051302 (2006).

[25] A Kudrolli, Size separation in vibrated granu-
lar material, Rep. Prog. Phys. 67, 209 (2004).

[26] J Duran, J Rajchenbach, E Clément, Arching
effect model for particle size segregation, Phys.
Rev. Lett. 70, 2431 (1993).

[27] J Duran, T Mazozi, E Clément, J Rajchen-
bach, Size segregation in a two-dimensional
sandpile: Convection and arching effects,
Phys.Rev. E 50, 5138 (1994).

[28] R O Uñac, A M Vidales, L A Pugnaloni, Sim-
ple model for wet granular beds subjected to
tapping, Gran. Matt. 11, 371 (2009).

[29] P Ribière, P Richard, P Philippe, D Bideau, R
Delannay, On the existence of stationary states
during granular compaction, Eur. Phys. J. E
22, 249 (2007).

070007-11