Journal Pre-proof

Shear-induced particle segregation in binary mixtures: Verification of a

percolation theory

Monica Tirapelle, Silvia Volpato, Andrea C. Santomaso

PII: S1674-2001(21)00037-7
DOI: https://doi.org/10.1016/j.partic.2021.01.005
Reference: PARTIC 1425

To appear in: Particuology

Received Date: 3 November 2020

Revised Date: 21 December 2020
Accepted Date: 5 January 2021

Please cite this article as: Monica Tirapelle, Silvia Volpato, Andrea C. Santomaso,
Shear-induced particle segregation in binary mixtures: Verification of a percolation theory,
<![CDATA[Particuology]]> (2021), doi: https://doi.org/10.1016/j.partic.2021.01.005

This is a PDF file of an article that has undergone enhancements after acceptance, such as
the addition of a cover page and metadata, and formatting for readability, but it is not yet the
definitive version of record. This version will undergo additional copyediting, typesetting and
review before it is published in its final form, but we are providing this version to give early
visibility of the article. Please note that, during the production process, errors may be
discovered which could affect the content, and all legal disclaimers that apply to the journal
pertain.

© 2020 Published by Elsevier.

 Shear-induced particle segregation in binary mixtures:
Verification of a percolation theory

Monica Tirapelle∗, Silvia Volpato, Andrea C. Santomaso∗∗

APTLab - Advanced Particle Technology Laboratory, Department of Industrial Engineering,
University of Padova, Via Marzolo 9, 35131, Padova, Italy

of
Abstract

ro
Granular materials composed of different-sized grains may experience undesired
segregation. Segregation is detrimental for a lot of industries because it leads to
an increase in production costs and wastes. For these reasons, the segregation

-p
phenomena have been intensively studied in the last decades, and a lot of model
have been provided by many researchers. However, these models are mainly
re
based on empirical relations rather than physical considerations. This paper
aims to confirm the main assumptions made by Volpato et al. (2020) in their
percolation theory by means of DEM simulations. The simulated geometry is a
lP

tilting shear box filled with few tracer particles in a bed of coarser sized grains,
and simulations are performed for a range of tilting frequencies and size ratios.
The results provide meaningful insight on the mathematical model parameters
na

and allow us to say that the percolation theory relies on physically consistent
assumptions.
Keywords: Discrete Element Method, Shear-induced percolation, Segregation,
ur

Parameter sensitivity analysis, Percolation model

Jo

∗ Corresponding author.

Email address : monica.tirapelle@phd.unipd.it

∗∗ Principal corresponding author.

Email address : andrea.santomaso@unipd.it

1
 1. Introduction

Granular materials may experience unexpected segregation owing to differ-

ences in size, shape, density and surface properties (Gray, 2018). Segregation
is typically described as a problem that plagues several process industries, such
as food, pharmaceutical, polymer but also mining, agricultural and chemical
industries (Hogg, 2009; Standish, 1985; Xiao et al., 2019). Segregation may

of
indeed determine hot spot and selectivity problems in reactors, may affect the
performance of a catalytic packed bed, may degrade the quality of a final prod-
uct (Bridgwater et al., 1985; Gray, 2018; Scott and Bridgwater, 1975). These

ro
are just a few examples of the many adverse situations causing an increase in
production costs and wastes.

-p
In the literature, at least 13 different mechanisms of segregation can be found
including percolation, kinetic sieving and squeeze expulsion. These mechanisms
however are often classified with some level of confusion. Also, the distinction
re
between causes and effects of segregation is sometimes unclear. To make an
example, the trajectory segregation mechanism should not be considered as a
lP

mechanism, but rather the result of differential drags and inertial forces between
different sized particles. Also, percolation is not really a basic segregation mech-
anism but the result of steric, frictional and inertial effects. In practice, it is
the combination of few basic mechanisms (steric, inertial, frictional, drag) that
na

give rise to very rich segregation phenomena that are, sometimes, improperly
referred to as segregation mechanisms. Aware of these ambiguities, in this paper
we will study the segregation phenomenon commonly known as percolation.
ur

Percolation was first suggested by Rosato et al. (1986) as responsible for

small particles to filter through the gaps in a matrix of larger grains. Percolation
Jo

occurs spontaneously if the particle size ratio, df /Dc , is smaller than 0.155
(Bridgwater and Ingram, 1971). For higher size ratios, it is commonly referred to
as shear-induced percolation or kinetic sieving and it is due to shear or vibrations
(Hogg, 2009; Khola and Wassgren, 2016). Due to the strain applied across the
failure zone, the larger particles will yield a space into which a smaller particle

2
 can move downward (Bridgwater, 1994; Scott and Bridgwater, 1975). Opposed
to kinetic sieving, there is squeeze expulsion. This phenomenon describes the
squeeze of a particle into an adjacent layer, also in the upward direction, and
affects all the particles, regardless of their properties. Combined with sieving,
squeeze expulsion determines the net percolation velocity of each species (Savage
and Lun, 1988).
A lot of models for describing segregation phenomena in different appara-

of
tuses were provided by many researchers in the last decades. However, because
they are mainly based on empirical relations rather than detailed physical con-

ro
siderations, a fundamental understanding of how grains percolate and segregate
is still lacking.
In this work, our main aim is to confirm the key assumptions made by Vol-

-p
pato et al. (2020) in their percolation theory by means of Discrete Element
Method (DEM) simulations. The simulated geometry is a shear box similar to
re
that used by Volpato et al. (2020), with free surface on the top and two moving
sidewalls. The latter are rotated backwards and forwards up to an inclination
angle of 45◦ and in such a way to generate a shear flow with a linear velocity pro-
lP

file (and hence a constant shear rate). The simulation domain has non periodic
boundary conditions and the bed of grains is subjected to constant strain rates.
DEM simulations are based on a non-linear spring-dashpot contact model: the
na

Hertz-Mindlin model with Constant Directional Torque (CDT) rolling friction.

This contact model requires the material properties (e.g. Young’s modulus,
friction coefficients...) to be defined. Since we could not directly obtain all the
ur

input parameters from physical experiments, we verified the robustness of our

parameters (and hence of our DEM simulations) through a parameter sensitiv-
ity analysis. The sliding friction coefficient is the only parameter significantly
Jo

influencing the results, and therefore had to be tuned accurately. Once tuned,
we used DEM simulations to provide meaningful insight on the mathematical
model parameters for the percolation theory (see the conceptual scheme de-
picted in Fig. 1). We found that the percolation theory in Volpato et al. (2020)
relies on physically consistent assumptions.

3
 of
ro
Fig. 1: Framework of the work. The light grey boxes concern a previously published work of
Volpato et al. (2020) whereas the dark ones are discussed in this paper.

-p
The paper is organized as follows. It starts by introducing the numeri-
cal method and describing the mathematical percolation model developed by
re
Volpato et al. (2020). Next, the percolation velocities predicted by the DEM
simulations are compared with experimental findings for parameter tuning. The
lP

validity of DEM simulations is strengthened by a parameter sensitivity analysis.

To conclude, the assumption of isotropic bed on which the mathematical model
is based has been checked numerically by characterising the volumetric porosity
by means of the Voronoi tessellation.
na

2. Numerical method
ur

2.1. Simulation method

One possibility to obtain information about the behaviour of granular media

is to simulate them with Discrete Element Method (Luding, 2008). The method
Jo

is simple, flexible, general and allows the forces acting on each grain to be
calculated by the integration of the Newton’s second law (Cundall and Strack,
1979).
If we consider two particles labelled i and j, the governing equations for
their translational and rotational motion and for their moment of inertia are

4
 respectively given by:
n
dvi X
mi = Fij + Fext,i , (1)
dt
j=1,j6=i

n
dωi X
Ii = (Ri × FijT ) + τij , (2)
dt
j=1,j6=i

where mi , Ri , Ii , vi and ωi are, in order, the mass, the radius, the moment of
inertia, the linear and the angular velocity of the i-th particle. τij is instead a

of
torque term that accounts for the effect of rolling friction during contact (Remy
et al., 2010). Notice that the contact force between the two grains, Fij , is

ro
usually decomposed into its normal and tangential components: FijN and FijT
whereas the external force, Fext,i , is most often represented by the force of
gravity: Fext,i = mi gi (Radjaı̈ and Dubois, 2011; Remy et al., 2010; Zhu et al.,
2007).

-p
To run the simulations, we employed the open source DEM particle simu-
®
re
lation software LIGGGHTS -PUBLIC (Kloss et al., 2012). As contact force
model, we used the non-linear spring-dashpot model developed by Hertz and
Mindlin (Johnson, 1987; Mindlin, 1949). So, if two spherical particles are in
lP

contact, namely their distance r is less than their contact distance d = Ri + Rj ,

the normal and tangential contact forces are calculated as:

FijN = kn δnij − γn vnij ,

na

(3)

FijT = kt δtij − γt vtij , (4)

where kn and kt are spring elastic constants, γn and γt are viscous damping
ur

constants, δnij and δtij are the normal and tangential displacements, and vnij
and vtij are relative velocities (LIGGGHTS(R)-PUBLIC website, n.d.). Fur-
thermore, the tangential force is governed by the Coulomb’s condition:
Jo

|FijT | ≤ µFijN , (5)

with µ being the friction coefficient. With this constraint, when the Coulomb
criterion is not satisfied, particle sliding occurs. The kn , kt , γn , γt coeffi-
cients are calculated from the material properties, so we had to set Young’s

5
 modulus, Poisson ratio, sliding friction coefficient and coefficient of restitution
(LIGGGHTS(R)-PUBLIC website, n.d.).
To account for gear-like rotation of two particles in contact, the contribution
of a rolling velocity is added in the contact model (Ai et al., 2011; Chialvo
et al., 2012; Radjaı̈ and Dubois, 2011). As rolling friction model we used the
Constant Directional Torque, CDT, in which the rolling friction is represented
by a constant torque, Mr , whose direction is the relative rotation between the

of
two bodies (LIGGGHTS(R)-PUBLIC website, n.d.) and reads:

ωi − ωj
Mr = − µr Rr Fn . (6)

ro
|ωi − ωj |

In Eq. 6, ωi and ωj are angular velocities and Rr is the rolling radius defined
as Rr ≡ ri rj /(ri + rj ) (Ai et al., 2011). Despite the CDT model can generate

-p
non stopping-oscillating torque in pseudo-static systems, it has been used here
because of its simple formulation and the low required computational effort.
re
Note that, in addition to the four parameters listed above, we had to define also
the coefficient of rolling friction.
lP

2.2. System geometry

The geometry used for generating and testing shear-induced percolation in

the granular assembly is a shear box similar to that used by Volpato and co-
na

workers in their experiments (Volpato et al., 2020). It consists of a transparent

bumpy glass bottom, two wooden fixed sidewalls and two wooden tilting side-
walls. This set-up bears several similarities with the shear box used by Bridg-
ur

water (1994) but unlike this one, the top surface is open and not confined, so
that the vertical stresses in the sheared bulk of grains are only due to the weight
of the grains and varies with depth (Royer and Chaikin, 2015). The box is 0.2
Jo

m deep and 0.2 m high, the distance between the tilting walls is 0.1 m wide. In
this system, gravity acts in the negative Z direction, whereas width and depth
are oriented respectively along the X and Y axes. A schematic representation
of the system is reported in Fig. 2.

6
 fig2.png

of
ro
-p
Fig. 2: A schematic representation of the simulated geometry in (a) its initial position and
re
(b) at the maximum inclination angle.

The tilting sidewalls remain parallel to each other and can rotate backward
lP

and forward until the wanted inclination angle of 45◦ (αmax ) in order to produce
a strain within the bulk of grains. To impose a constant shear rate, (γ̇ ≡
v(h)/h = const), the angular velocity, ω, imposed to the rotating walls reads:
na

arctan(γ̇ · t)
ω(t) = , (7)
t

with 0 < t < t(45◦ ). It is worth noticing that, despite the shear rate was held
ur

constant throughout most of the forward and backward movement of the cell,
transients necessarily occur during the direction reversal (Khola and Wassgren,
2016).
Jo

The tilting box, unlike other types of shearing devices, has the disadvantage
of a periodic discontinuous motion. A continuous and indefinite shearing action
can be provided, for instance, in an annular shear cell (Artoni et al., 2018; May
et al., 2010; Savage and Sayed, 1984; Wildman et al., 2008). In the latter case
however, the velocity field is strongly dependent on the distance from the bottom

7
 because an inhomogeneous shear rate profile would develop in the granular
assembly. On the contrary, the tilting box can provide a homogeneous (constant)
shear in the granular bed, with the sole exception of the upper portion of the
bed where some recirculation cells develop (as discussed in detail in Sect. § 2.3).
Our simulated shear box is different from the other shear boxes simulated
and reported in the recent literature. For instance, unlike van der Vaart et al.
(2015), we applied constant shear rate and, differently from Khola and Wassgren

of
(2016), walls were simulated as frictional and boundaries as non-periodic. These
choices were done in order to mimic the experimental shear box.

ro
2.3. Particle properties and input parameters

A full factorial design of simulations was performed for the following size

-p
ratios: d∗ = df /Dc =0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55 and for the following
three different tilting frequencies: 17, 26 and 35 rpm, which correspond respec-
re
tively to shear rates of 0.57, 0.87 and 1.17 s−1 . The 21 simulated combinations
are reported in Table 2.
The cell was filled with spherical glass beads assumed to be monodisperse
lP

and cohesionless, up to a bed height of 16 cm. The bed of coarse particles

was generated within the cell and let settle by gravity. In order to prevent the
development of a regular hexagonal close-packing structure close to the bottom
na

(Scott and Bridgwater, 1975), the first layer of particles presented a polydisperse
size distribution (Gaussian distribution with µ = 6 mm and σ = 0.15 · µ. The
standard deviation is big enough to avoid regular packing, but it likewise ensures
ur

that the particles in the bed are always larger than the tracer particles); whereas
all the other grains in the bulk had diameter equals to 6 mm. The insertion
of fine percolating particles was successively done by randomly replacing some
Jo

of the coarse grains placed at a constant distance from the bottom (10 · Dc ),
spaced each other in order to ensure very diluted conditions (i.e. few fines in a
bed of coarse particles), and far enough (at least 2 Dc far) from the sidewalls
in order to minimize any wall influence. The fine percolating particles were
inserted below the free surface in order to provide a sufficiently thick layer of

8
 coarse particles above them. This layer prevents the fine particles to be trapped
into the recirculation loops that typically develop in the upper part of the bed
below the free surface. Indeed Fig. 3 shows that significant deviations from the
imposed vx component appear only above 12 cm bed height (20 · Dc ), regardless
of the strain applied. The same investigation was done also for a 12 cm bed
height (20 · Dc ) and even in that case, only the grains in the lower 75% of bed
experienced a shear in horizontal direction. The fine percolating particles were

of
positioned with the same care also in the experimental study (Volpato et al.,
2020).

ro
Fig. 3: The horizontal velocity component (vx ) of particles considering a time lapse of 1
second during which the bed is sheared, in the case of 17, 26 and 35 rpm tilting frequencies.

-p
The dashed orange lines represent the theoretical vx whereas the black ones denote the limit
below which undisturbed shear is present.
re
Bulk particles and percolating particles had the same density of 2540 kg/m3 .
The total amount of grains in each simulation was around 16000. Once the
material was loaded and the steady state attained, the sidewalls started to tilt.
lP

The material properties and the simulation parameters of the standard case
are summarized in Table 1. All properties are typical of glass beads, except for
the elastic Young’s modulus. Since it has been shown that the use of a relatively
na

small Young’s modulus does not result in a significant error in structural anal-
ysis (Zhou et al., 2004), E was set of the order of magnitude of MPa instead of
GPa. This reduces the computational time without significantly affecting flow
ur

patterns, velocity profiles and shear stresses (LIGGGHTS(R)-PUBLIC website,

n.d.; Remy et al., 2009). The wall-particle friction was instead imposed equal to
0.18 regardless of the type of materials in contact because µwp was experimen-
Jo

tally measured equal to 0.17 ± 0.01 for wood and glass, and equal to 0.19 ± 0.01
for glass and glass. However, because the input parameters play an important
role in generating accurate results (Zhu et al., 2007), the parameter tuning was
successively checked by implementing a sensitivity analysis.
The numerical time step, which should be smaller than a critical value for

9
 Variable Symbol Value
Particle density (kg/m3 ) ρs 2540
Young’s modulus (MPa) E 26
Poisson ratio σ 0.25
Sliding friction coefficient µs 0.18
Rolling friction coefficient µr 0.005
Wall-particle friction µwp 0.18
Restitution coefficient en 0.60

of
Coarsest particle diameter (m) Dc 0.006
Diameter ratio df /Dc 0.25-0.55
Number of particles NT > 16000

ro
Time step (s) ∆t < 1x10−5
Shear rate (s−1 ) γ̇ 0.567-1.167
Boundary conditions − fff

-p
Table 1: A summary of the DEM simulation parameters of the standard case.

avoiding, in each time step, the propagation of the disturbance farther than each
re
particle immediate neighbourhood (Zhu et al., 2007), was chosen to be at least
20% lower than the Rayleigh time. In particular ∆T was imposed 0.5 · 10−5 s
for d∗ = 0.25 − 0.30, and 1.0 · 10−5 s in all the other cases.
lP

3. Mathematical model
na

This section briefly describes the mathematical model used to fit the data
obtained from numerical simulations. The model was proposed by Volpato et al.
(2020) and allows the percolation velocity of a fine isolated spherical particle
ur

in a sheared bed of coarser monomodal grains to be determined. The model

considers only the contribution of kinetic sieving and disregards the squeeze
expulsion since it is supposed that squeeze expulsion, in the absence of signifi-
Jo

cant density difference between particles, contributes equally in the upward and
downward directions, affecting minimally the overall percolation velocity. It is
indeed assumed that squeeze expulsion is not size preferential and that there is
no inherent preferential direction for the layer transfer (Savage and Lun, 1988).
According to the model, which arises from statistical considerations under the

10
 hypothesis of isotropic bed dilation, the dimensionless percolation velocity is:
A · (1 − d∗ ) · Pf
vp∗ = √ , (8)
Pf · γ̇ ∗ 1 − d∗ + (1 − Pf )
with the three dimensionless variables vp∗ , γ̇ ∗ and d∗ being respectively equal to:
q
v d
vp∗ = γ̇Dpc , γ̇ ∗ = γ̇ Dgc , d∗ = Dfc .

Pf represents the probability for the fine particle to find an aperture large

of
enough to fall-in under the action of gravity and it was demonstrated to be:
 
3 df (1 − ε)
Pf = exp − , (9)
k Dc ε

ro
where ε is the bed porosity. As it can be seen, the model has two parameters: A
and k. The former is a proportionality parameter expected to be dependent on

-p
the material properties, whereas the latter depends on the geometrical properties
of the particle bed and it is equal to 2/3 in the case of isotropic bed under static
conditions (Snabre and Mills, 2000; Volpato et al., 2020). Detailed information
re
on the percolation model can be found elsewhere (Volpato et al., 2020).
lP

4. Results and discussion

4.1. Parameter tuning

In order to have realistic simulation results, the parameters need to be cal-
na

ibrated. However, calibration has a problem of ambiguity of the parameters’

selection (Roessler et al., 2019). Furthermore, it is difficult to obtain all the
parameters by physical experiments. Therefore, we performed a parameter sen-
ur

sitivity analysis to evaluate the robustness of the parameters we set. This study
was done by systematically changing the mechanical properties of glass beads
and evaluating their own impact on the percolation velocity, considering differ-
Jo

ent size ratios and 26 rpm as tilting frequency. Once evaluated the percolation
velocities, a response surface regression was made for each parameter. Note that
the effect of the mutual interaction of the parameters has not been studied.
The coefficient of rolling friction is usually added to mimic the behaviour of
not perfectly spherical particles. Here, we varied the rolling friction in the range

11
 Run Rpm df /Dc

of
1 17 0.25
2 17 0.30
3 17 0.35

ro
4 17 0.40
5 17 0.45
6 17 0.50

-p
7 17 0.55
8 26 0.25
9 26 0.30
10 26 0.35
re
11 26 0.40
12 26 0.45
13 26 0.50
lP

14 26 0.55
15 35 0.25
16 35 0.30
17 35 0.35
na

18 35 0.40
19 35 0.45
20 35 0.50
21 35 0.55
ur

Table 2: The combination treatments of the full factorial design of experiments. The columns
are: ID number, tilting frequency in rpm and size ratio.
Jo

12
 (a) (b)

(c) (d)

Fig. 4: Dimensionless percolation velocity as a function of the size ratio in the case of 26
rpm tilting frequency and for different values of (a) rolling friction, (b) Young’s modulus, (c)
restitution coefficient and (d) Poisson ratio.

of
of 0.005-0.100 which approximately match the one of glass beads. As Fig. 4-a
shows, the rolling friction has no drastic influence on the percolation velocity.

ro
Neither Young’s modulus in the range of 1 · 107 -5 · 107 Pa, nor the coefficient
of restitution, which was varied in between 0.4 and 0.8, have significant impact
on the results (see Fig. 4-b,c). The same conclusion can be drawn also for

-p
variations in the Poisson ratio from 0.20 to 0.30. (Fig. 4-d). According to
these simulations, the shape and the mechanical properties of the material have
re
minimal influence on the percolating velocity. In each case indeed, it appears
from response surface regression with backwards elimination that the response
variable (namely the dimensionless percolation velocity) depends quadratically
lP

on the size ratio but it is independent on the parameter under consideration.

What strongly influence the granular behaviour of the system are indeed the
frictional properties. As Fig. 5 shows, higher frictional coefficients determine an
na

increase in the percolation velocity. In this case the response has a statistically
significant dependence on both sliding friction coefficient, size ratio and their
interaction:
ur

vp∗ = 8.63 − 34.18d∗ + 5.09µs + 35.33d∗2 − 10.33d∗ µs (10)

The standard deviation of the data points around the fitted values, which is
Jo

equal to 0.20 (in dimensionless units since vp∗ is dimensionless), and R2 = 96.5%
ensure the reliability of the model. Fig. 6 represents the corresponding contour
map. It is interesting to see that for a fixed d∗ the percolation velocity in-
creases with friction. Despite this result can be counter-intuitive at a first sight,
different explanations can be found in the available literature. For instance,

13
 according to Remy et al. (2009), this is due to an increase in the granular tem-
perature with associated an increase in diffusive mixing; for Jing et al. (2017),
increasing interparticle friction promotes the upward migration of large particles
and hence, the creation of voids that are easily filled by adjacent fines (Hogg,
2009). According to us, the effect of friction on the percolation process is prob-
ably related to changes in the packing structure of the sheared bed. In order
to verify this hypothesis, we performed dedicated simulations on the shear cell

of
measuring the bed porosity as a function of the sliding friction coefficient. Fig.
7 shows an increase of the bulk porosity with increasing the friction coefficient.

ro
The lower mobility of the coarse particles due to the increased friction at the
contact points determines a less efficient packing during the bed deposition and
also during the shearing action. According to Eq. 9, this implies larger proba-

-p
bility for a percolating particle to find a void large enough to percolate through
and hence, larger percolation velocity. Clearly, the sliding friction coefficient is
re
the only parameter that requires an accurate tuning. We set µs = 0.18. The
adoption of µs = 0.18 is justified both from direct experimental measurements
(µs = 0.19 ± 0.01) and because µs = 0.18 corresponds exactly to the estimated
lP

bulk porosity of ε = 0.40 (see Fig. 7).

The effect of changes in packing quality can also be observed in Fig. 8 that
reports the trends of the parameters A and k (present in Eq. 9) as a function
na

of bed porosity. k is not significantly affected by the average bed porosity and,
except for the lowest value, which corresponds to a very low interparticle friction
(Fig. 7), it tends to a constant value. This is also physically accurate because
ur

k is a geometrical parameter and hence expected to be more affected by some

shear-induced anisotropy rather than changes in the average bed porosity. On
the other hand, the parameter A increases linearly with porosity, suggesting
Jo

the existence of a relationship between A and the average distance between the
particles. However, this effect is already incorporated into the probability term
(Eq. 9) and hence, we here suggest a dependence of A on interparticle friction
mediated by the bulk porosity. This is also shown in Fig. 5 where we observe
a monotonic increase of A with the interparticle friction coefficient (i.e. with

14
 Fig. 5: Effect of changes in sliding friction coefficient on the percolation velocity. The tilting
frequency used is 26 rpm.

Fig. 6: Contour plot of dimensionless percolation velocity versus size ratio and friction coeffi-
cient.

a porosity increment). A larger distance between particles increases the falling

of
distance and therefore the velocity of the fine percolating particles trapped in
the cage made of coarse ones, increasing the overall percolation velocity.

ro
4.2. Percolation velocities

Simulations were performed for a range of size ratios and strain rates. The

-p
resulting combinations of the full factorial design of experiments are reported in
Table 2. The average percolation velocity, which is the response variable, was
re
calculated for each run as:
N
P
vp,j
j=1
vp = , (11)
N
lP

where vp,j is the percolation rate of the j-th particle to cross the cell and N is
the number of fine percolating grains. N was 32 for d∗ = 0.20, 0.25, 0.30, and
16 for all the other cases.
na

Fig. 9 shows the velocity profile as a function of size ratio for each tested
shear rate when using the parameters listed in Table 1 (remind that the relia-
bility of these parameters has been proved in §4.1). The error bars represent
ur

the standard deviation of 16 or 32 independent measurements, which becomes

larger for smaller percolating particles. As expected, the percolation velocity
significantly decreases for size ratios tending to 1 (at 5% significance level) in-
Jo

deed, when percolating particles are surrounded by particles with size closer to
its own size, the segregation flux diminishes. On the other hand, the smaller

Fig. 7: Porosity as a function of the coefficient of sliding friction.

15
 Fig. 8: The parameters A (on the left axis) and k (on the right axis) are plotted as a function
of the bed porosity.

Fig. 9: Percolation velocity as a function of the size ratio for three different frequencies of
tilting: 17, 26 and 35 rpm. In the inset, the dimensionless percolation velocities are reported.

the falling particles, the higher number of opportunities to find void bigger than

of
themselves and the longer the falling distance. The percolation velocity is larger
at high shear rates.

ro
In the inset of Fig. 9, the same data are represented in terms of dimensionless
percolation velocity. As found in Volpato et al. (2020), all data collapse on the

-p
same curve since the dimensionless percolation velocity is independent of the
shear rate.
In Fig. 10 the simulated dimensionless velocities are fitted with the mathe-
re
matical model presented in Section §3 by using the Least-squares minimization
(in detail: Trust Region Reflective minimization algorithm (Virtanen et al.,
lP

2020)). The optimal values found for the parameters, considering ε = 0.40, are:
A = 10.80 and k = 0.80. The latter suggests that, as the size ratio approaches
0.80, the percolation becomes negligible (the velocity profile converges to 0). It
should be noted that k = 0.80 is slightly larger than 2/3, typical value that
na

characterizes isotropic static beds. The root-mean-square deviation (RMSD)

between predicted and observed values is 0.133, showing a good quality of the
fitting.
ur

The numerical results reveal also a nice agreement with the previously re-
ported experimentally achievements (Volpato et al., 2020). In that case, A and
Jo

k were respectively equal to 9.70 and 0.71 and hence, consistent with what is
found here. The numerical model is therefore valid and can be used to accu-

Fig. 10: Fitting of the numerical data with the mathematical model proposed by Volpato
et al. (2020).

16
 Fig. 11: Distribution (on the left) and cumulative distribution (on the right) of percolation
velocities. The analysed particles are 176 in the case of 26 rpm tilting frequency and 0.45 size
ratio.

rately predict percolation in sheared systems with dilute fines concentration in

the range of the explored shear rates.

of
A typical percolation velocity distribution for 176 individual spherical parti-
cles is shown in Fig. 11. In this case the bed was sheared with a tilting frequency
of 26 rpm and d∗ = 0.45. The result is a left-skewed distribution meaning that

ro
more than half of the particles percolate faster than the mean percolation rate.

-p
4.3. Packing structure characterisation

The mathematical model (9) shows that the packing structure plays an im-
portant role on the percolation velocity: the probability Pf is indeed a function
re
of both local porosity ε and bed stereo-geometry k. It is therefore important to
characterize the internal packing structure of the granular bed.
lP

The bed porosity was computed by using the tessellation method based on
the 3D generalization of the Voronoi diagrams. By definition, a Voronoi cell
around a particle is the region of space that is much closer to that particle than
to any other particle in the system (Burtseva and Werner, 2015; The Royal
na

Society of Chemistry, 2010). So, the Voronoi tessellation divides the space into
regular polyhedra with flat faces and straight edges and allows the local porosity
to be calculated as:
ur

VV oro − Vp
εl = (12)
VV oro
In Eq. 12, VV oro and Vp are respectively the volume of the tessellation
Jo

containing the i-th particle and the volume of the particle itself. The latter
was known, whereas VV oro was computed in LIGGGHTS -PUBLIC by using ®
Voro++, an open source software library (Rycroft, 2009).
In Fig. 12 the density distribution and the cumulative distribution of the
local porosities are reported both in static conditions and in the case of sheared

17
 bed. The packing condition varies after the shear is applied: the bed porosity
slightly decreases from an initial median value of 0.42 to 0.40 indicating a con-
traction (settling) of the bed under shear with respect to the initial packing.
Note that the medians are more representative than the means because less bi-
ased by the outliers. The second peak of the distributions is indeed due to the
over-relaxation of the bed close to the walls while the peak at ε=1 refers to the
particles at the free surface of the bed. Both these peaks are not representative

of
of the bulk. The porosity ε = 0.40, which was used in the mathematical model
for fitting the simulation data (Fig. 10), is therefore justified and in agreement

ro
with the experimental findings of Volpato et al. (2020).
The porosity estimated from the Voronoi tessellation is a volumetric porosity
that does not provide any information about the isotropic or anisotropic nature

-p
of the packing of spheres (Grattoni and Dawe, 1995). Because an anisotropic
packing bed could affect the percolation process, and because the model (Vol-
re
pato et al., 2020) rests on the hypothesis of perfect isotropy of the bed, the
structure of our bed of coarse glass beads was studied in detail. To establish
if there was anisotropy, a cubic region of edge length 8 cm within the powder
lP

bed was considered. For visual reference, Fig. 13 shows the selected region and
examples of 2-D cross sections orthogonal to the X, Y and Z directions. The
black colour indicates the solid particles whereas the white colour the empty
na

spaces. The porosity characteristics of XY, XZ, and YZ planes is calculated as

the ratio between the empty spaces (white pixel) and the area of the whole cross
section. In order to calculate a porosity distribution for each of the three tested
ur

directions, porosity was estimated on a hundred of slices along any direction. In

Fig. 14 the respective boxplots are shown. Unlike along x and y directions, the
porosity along z spreads over a wider range and some outliers appear. However,
Jo

as detailed below, this difference is not statistically significant.

The assumption of isotropic packing structure on which Eq. 9 relies is justify
below by means of statistics. Firstly, the Bartlett’s test of equality of variances
was performed at the 0.05 significance level. Because the p-value is less than
the significance level, the hypothesis of equal variances has to be rejected: a

18
 Fig. 12: Distribution (on the left) and cumulative distribution (on the right) of the bed
porosity in static (blue bars) and sheared condition (orange bars).

statistically significant difference exists between the variances of at least two

independent sets of our normally distributed continuous data. Then, for testing
the equality of the porosity in the three directions, a one-way ANOVA F-test

of
statistics scheme was applied assuming no equal variances. Because the p-
value=0.721 is greater than the significance level of 5%, there is no statistically

ro
significant differences between the means. We can therefore conclude that the
mean porosity is the same in all three directions and hence, the assumption of
isotropic bed under shear conditions originally made in Volpato et al. (2020) is
justified.

-p
re
5. Conclusions

In this work, DEM simulations of a shear box were used to study shear
lP

induced percolation. In particular, we focused on the segregation velocity of

few small particles in a bed of coarser sized grains. The simulations spanned a
range of size ratios and a range of shear rates that are characteristic of the quasi-
static regime of flow (Tardos et al., 2003), namely until dimensionless shear rates
na

γ̇ ∗ of 0.15–0.20. As already reported in literature (Bridgwater, 1994; Khola and

Wassgren, 2016; Volpato et al., 2020), it was found that the percolation speed
decreases exponentially when the particle size ratio increases; furthermore the
ur

dimensionless percolation velocity is independent of the shear rate, at least for

shear rates in the quasi-static regime of flow.
Jo

The reliability of our DEM simulations is proved by a parameter sensitivity

analysis. Parameters such as Young’s modulus and Poisson ratio do not in-
fluence the results, unlike the sliding friction coefficient. However, the sliding
friction coefficient is tuned correctly as evidenced by the experiments.
After the parameter tuning, DEM simulations were used to verify the hy-

19
 of
ro
-p
re
lP
na
ur
Jo

Fig. 13: Orthogonal slices for 3D visualisation of the pore structure.

20
 Fig. 14: The boxplots of the bed porosity along the three directions: x, y and z.

pothesis of an isotropic bed dilation on which the mathematical percolation

model of Volpato et al. (2020) is based. In sheared conditions, the pack of
grains results to be isotropic with volumetric median porosity of 0.40 and fric-
tion coefficient equals to µs = 0.18. Also, the parameters A and k are in close

of
agreement with the percolation theory and rely on physical considerations: k is
related to some shear-induced anisotropy, whereas A depends on interparticle
friction mediated by bulk porosity.

ro
We can conclude that the percolation theory developed by Volpato et al.
(2020) relies on physically consistent assumptions and that the mathematical

-p
model parameters are meaningful. A future natural extension for the percolation
theory would be to consider the effect of fine particles concentration, and to
introduce polydispersity and particle density differences. The validation of such
re
a model would be difficult by means of physical experimentation. However, since
we have already tuned the parameters, we could rely on DEM simulations.
lP

References

Ai, J., Chen, J.F., Rotter, J.M., Ooi, J.Y., 2011. Assessment of rolling resistance
na

models in discrete element simulations. Powder Technology 206, 269–282.

Artoni, R., Soligo, A., Paul, J.M., Richard, P., 2018. Shear localization and
wall friction in confined dense granular flows. Journal of Fluid Mechanics
ur

849, 395–418.

Bridgwater, J., 1994. Mixing and segregation mechanisms in particle flow, in:
Jo

Granular matter. Springer, pp. 161–193.

Bridgwater, J., Cook, H.H., Drahun, J.A., 1985. Strain induced percolation, in:
Instn. Chem. Engrs. Symposium Series, pp. 171–191.

21
 Bridgwater, J., Ingram, N.D., 1971. Rate of spontaneous inter-particle percola-
tion. Trans. Inst. Chem. Engineer 49, 163–169.

Burtseva, L., Werner, F., 2015. Modeling of spherical particle packing structures
using mathematical tessellation. Univ., Fak. für Mathematik.

Chialvo, S., Sun, J., Sundaresan, S., 2012. Bridging the rheology of granular
flows in three regimes. Physical review E 85, 021305.

of
Cundall, P.A., Strack, O.D., 1979. A discrete numerical model for granular
assemblies. Geotechnique 29, 47–65. doi:10.1680/geot.1979.29.1.47.

ro
Grattoni, C.A., Dawe, R.A., 1995. Anisotropy in pore structure of porous media.
Powder technology 85, 143–151.

-p
Gray, J.M.N.T., 2018. Particle segregation in dense granular flows. Annual
Review of Fluid Mechanics 50, 407–433.
re
Hogg, R., 2009. Mixing and segregation in powders: evaluation, mechanisms
and processes. KONA Powder and Particle Journal 27, 3–17.
lP

Jing, L., Kwok, C.Y., Leung, Y., 2017. Micromechanical origin of particle size
segregation. Physical review letters 118, 118001.
na

Johnson, K.L., 1987. Contact mechanics. Cambridge university press.

Khola, N., Wassgren, C., 2016. Correlations for shear-induced percolation seg-
regation in granular shear flows. Powder technology 288, 441–452.
ur

Kloss, C., Goniva, C., Hager, A., Amberger, S., Pirker, S., 2012. Models, algo-
rithms and validation for opensource dem and cfd–dem. Progress in Compu-
Jo

tational Fluid Dynamics, an International Journal 12, 140–152.

LIGGGHTS(R)-PUBLIC website, n.d. https://www.cfdem.com.

Luding, S., 2008. Introduction to discrete element methods: basic of contact

force models and how to perform the micro-macro transition to continuum

22
 theory. European Journal of Environmental and Civil Engineering 12, 785–
826.

May, L.B., Golick, L.A., Phillips, K.C., Shearer, M., Daniels, K.E., 2010. Shear-
driven size segregation of granular materials: Modeling and experiment. Phys-
ical Review E 81, 051301.

Mindlin, R.D., 1949. Compliance of elastic bodies in contact. Journal of Applied

of
Mechanics 16, 259–268.

Radjaı̈, F., Dubois, F., 2011. Discrete-element modeling of granular materials.

ro
Wiley-Iste.

Remy, B., Canty, T.M., Khinast, J.G., Glasser, B.J., 2010. Experiments and

-p
simulations of cohesionless particles with varying roughness in a bladed mixer.
Chemical Engineering Science 65, 4557–4571.
re
Remy, B., Khinast, J.G., Glasser, B.J., 2009. Discrete element simulation of
free flowing grains in a four-bladed mixer. AIChE Journal 55, 2035–2048.
lP

Roessler, T., Richter, C., Katterfeld, A., Will, F., 2019. Development of a
standard calibration procedure for the DEM parameters of cohesionless bulk
materials – part I: Solving the problem of ambiguous parameter combina-
na

tions. Powder Technology 343, 803–812. URL: https://doi.org/10.1016/

j.powtec.2018.11.034, doi:10.1016/j.powtec.2018.11.034.

Rosato, A., Prinz, F., Standburg, K.J., Swendsen, R., 1986. Monte carlo simu-
ur

lation of particulate matter segregation. Powder Technology 49, 59–69.

Royer, J.R., Chaikin, P.M., 2015. Precisely cyclic sand: Self-organization of

Jo

periodically sheared frictional grains. Proceedings of the National Academy

of Sciences 112, 49–53.

Rycroft, C., 2009. Voro++: A three-dimensional Voronoi cell library in C++.

Technical Report. Lawrence Berkeley National Lab.(LBNL), Berkeley, CA
(United States).

23
 Savage, S.B., Lun, C.K.K., 1988. Particle size segregation in inclined chute flow
of dry cohesionless granular solids. Journal of Fluid Mechanics 189, 311–335.

Savage, S.B., Sayed, M., 1984. Stresses developed by dry cohesionless granular
materials sheared in an annular shear cell. Journal of Fluid Mechanics 142,
391–430.

Scott, A.M., Bridgwater, J., 1975. Interparticle percolation: a fundamental

of
solids mixing mechanism. Industrial & Engineering Chemistry Fundamentals
14, 22–27.

ro
Snabre, P., Mills, P., 2000. Settling and fluidization of non Brownian hard
spheres in a viscous liquid. The European Physical Journal E 1, 105–114.

-p
doi:10.1007/PL00014590.

Standish, N., 1985. Studies of size segregation in filling and emptying a hopper.
re
Powder Technology 45, 43–56.

Tardos, G.I., McNamara, S., Talu, I., 2003. Slow and intermediate flow of
lP

a frictional bulk powder in the couette geometry. Powder Technology 131,

23–39.

The Royal Society of Chemistry, 2010. Using a voronoi tessellation to measure

na

properties of a polydisperse sphere packing in a microphase separated system

.

van der Vaart, K., Gajjar, P., Epely-Chauvin, G., Andreini, N., Gray, J.M.N.T.,
ur

Ancey, C., 2015. Underlying asymmetry within particle size segregation.

Physical review letters 114, 238001.
Jo

Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Courna-
peau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al., 2020.
SciPy 1.0: Fundamental Algorithms for ScientificComputing in Python. Na-
ture Methods 17, 261–272.

24
 Volpato, S., Tirapelle, M., Santomaso, A.C., 2020. Modeling and experimental
investigation of shear-induced particle percolation in diluted binary mixtures.
Physical Review E 102, 012902.

Wildman, R.D., Martin, T.W., Huntley, J.M., Jenkins, J.T., Viswanathan, H.,
Fen, X., Parker, D.J., 2008. Experimental investigation and kinetic-theory-
based model of a rapid granular shear flow. Journal of Fluid Mechanics 602,

of
63–79.

Xiao, H., Fan, Y., Jacob, K.V., Umbanhowar, P.B., Kodam, M., Koch, J.F.,

ro
Lueptow, R.M., 2019. Continuum modeling of granular segregation during
hopper discharge. Chemical Engineering Science 193, 188–204.

-p
Zhou, Y.C., Yu, A.B., Stewart, R.L., Bridgwater, J., 2004. Microdynamic anal-
ysis of the particle flow in a cylindrical bladed mixer. Chemical Engineering
Science 59, 1343–1364. doi:10.1016/j.ces.2003.12.023.
re
Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2007. Discrete particle simula-
tion of particulate systems: theoretical developments. Chemical Engineering
lP

Science 62, 3378–3396.

na
ur
Jo

25
Graphical Abstract (for review) Click here to access/download;Graphical Abstract (for
review);graphical_abstract_shear_box.eps

of
ro
-p
re
Experiments

lP quali cation validation

Mathematical Discrete Element
percolation model Method simulations
na
Hyp: isotropic bed dilation Parameters' adoption

veri cation of

model assumption
ur

Paramater
sensitivity analysis
Jo
Highlights

Highlights:

 Shear-induced percolation is studied with DEM.

 Parameter sensitivity analysis highlights the role of sliding friction on system response.

 The assumptions on which the mathematical percolation model relies are verified.

 The parameters of the mathematical model result to be physically consistent.

 Statically significant evidence of bed isotropy in sheared conditions is reported.

of
ro
-p
re
lP
na
ur
Jo
Figure 1 Click here to access/download;Figure;fig 01.png

of
ro
-p
re
lP
na
ur
Jo
Figure 2 Click here to access/download;Figure;fig 02.jpg

of
ro
-p
re
lP
na
ur
Jo
Figure 3 Click here to access/download;Figure;fig 03.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 4 Click here to access/download;Figure;fig 04.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 5 Click here to access/download;Figure;fig 05.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 6 Click here to access/download;Figure;fig 06.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 7 Click here to access/download;Figure;fig 07.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 8 Click here to access/download;Figure;fig 08.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 9 Click here to access/download;Figure;fig 09.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 10 Click here to access/download;Figure;fig 10.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 11 Click here to access/download;Figure;fig 11.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 12 Click here to access/download;Figure;fig 12.tif

of
ro
-p
re
lP
na
ur
Jo
Figure 13 Click here to access/download;Figure;fig 13.png

of
ro
-p
re
lP
na
ur
Jo
Figure 14 Click here to access/download;Figure;fig 14.tif

of
ro
-p
re
lP
na
ur
Jo
Declaration of Interest Statement

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships
that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered
as potential competing interests:

of
ro
-p
re
lP
na
ur
Jo