CHINA PARTICUOLOGY Vol. 3, No.

6, 354-363, 2005

DISCRETE AND CONTINUUM MODELLING OF

GRANULAR FLOW
H. P. Zhu1, Y. H. Wu2 and A. B. Yu1,*
1
Centre for Simulation and Modelling of Particulate Systems, School of Materials Science and Engineering,
the University of New South Wales, Sydney, NSW 2052, Australia
2
Department of Mathematics and Statistics, Curtin University of Technology, Perth, WA 6845, Australia
*Author to whom correspondence should be addressed. E-mail: [Link]@[Link]

Abstract This paper analyses three popular methods simulating granular flow at different time and length scales:
discrete element method (DEM), averaging method and viscous, elastic-plastic continuum model. The theoretical models
of these methods and their applications to hopper flows are discussed. It is shown that DEM is an effective method to
study the fundamentals of granular flow at a particle or microscopic scale. By use of the continuum approach, granular
flow can also be described at a continuum or macroscopic scale. Macroscopic quantities such as velocity and stress can
be obtained by use of such computational method as FEM. However, this approach depends on the constitutive rela-
tionship of materials and ignores the effect of microscopic structure of granular flow. The combined approach of DEM and
averaging method can overcome this problem. The approach takes into account the discrete nature of granular materials
and does not require any global assumption and thus allows a better understanding of the fundamental mechanisms of
granular flow. However, it is difficult to adapt this approach to process modelling because of the limited number of parti-
cles which can be handled with the present computational capacity, and the difficulty in handling non-spherical particles.
Further work is needed to develop an appropriate approach to overcome these problems.
Keywords granular flow, hopper flow, discrete element method, average method, continuum approach

1. Introduction DEM-based simulation has been used to study the physics

of granular materials (e.g., Herrmann & Luding, 1998; Ki-
The dynamic behavior of granular materials is very shino, 2001; Zhu et al., 2004), and applied to investigate
complicated. As with solids, they can withstand deforma- various engineering problems such as mixing (Wightman
tion and form heaps; as with liquids, they can flow; as with et al., 1998; Stewart et al., 2001), grinding/granulation
gases, they exhibit compressibility. These features give (Yamane et al., 1998; Yang et al, 2003) and hopper flow
rise to another state of matter that is poorly understood (de (Ristow, 1992; Langston et al., 1995a; Zhu & Yu, 2004).
Gennes, 1999). Corresponding to the fluid- and solid-like By use of a proper averaging procedure, the discrete
modes, different regimes have been identified in the past: system considered above can be transferred into a corre-
quasi-static regime, rapid flow regime and a transitional sponding continuum system. Extensive research has been
regime that lies in between. To describe the flow behavior carried out to develop such averaging methods. The
in these regimes, different approaches at different time and methods thus far proposed can be divided into three
length scales have been developed in the past. Generally classes: volume average (Drescher & de Josselin de Jong,
speaking, they can be categorize into two length scales: 1972; Rothenburg & Selvadurai, 1981), time-volume av-
microscopic and macroscopic. At a microscopic or particle erage (Walton & Braun, 1986; Zhang & Campbell, 1992)
scale, granular material is a discrete system whose and weighted time-volume average (Babic, 1997; Zhu & Yu,
physical properties are discontinuous with respect to posi- 2002). Based on a proper averaging method, the macro-
tion and time. On the other hand, at a macroscopic or bulk scopic quantities such as density, velocity and stress can
scale, granular material is a continuum system whose be obtained in terms of the microscopic quantities such as
physical properties are continuous. Based on such con- velocities of particles, and interaction forces and torques
siderations, different approaches have been proposed to between particles. The combined approach of DEM and
describe granular material at the two distinct scales. averaging method has been applied to the macro-dy-
An important simulation approach which can describe namical analysis of granular flows (e.g., Langston et al.,
granular flow at a particle scale is the so-called discrete 1995a; Potapov & Campbell, 1996; Zhu & Yu, 2005b), and
element method (DEM) originally developed by Cundall has been used to study the intrinsic characteristics of
and Strack (1979). The method considers a finite number granular materials such as the constitutive relationship
of discrete interacting particles, and every particle in the under various flow conditions (Oda & Iwashita, 2000;
considered system is described by Newton’s equations of Alonso-Marroquin & Herrmann, 2005).
motion related to translational and rotational motions. The Granular material can also be modeled by continuum
motion of every particle is traced and the gradients of approach on a macroscopic scale. In the continuum ap-
translation and rotation of a particle are determined in proach, the macroscopic behavior of granular flow is de-
terms of the forces and the torques exerted on it. scribed by the balance equations facilitated with constitu-
 Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow 355

tive relations and boundary conditions. In the past, two force is described by the dynamic friction, usually, given by
continuum models developed based on the plasticity the- the Coulomb friction model. The contact between two par-
ory and kinetic theory of molecular dynamics have exten- ticles is not a single point but is a finite area due to the
sively been used to study the dynamic behavior of granular deformation of both particles. The interparticle forces act
materials (Lun et al., 1984; Campbell, 1990; Nedderman, over the contact region between particles rather than the
1992). They have been illustrated to be applicable to mass centre of the particles, and they will generate a
quasi-static and rapid flow regimes, respectively. Recently, torque. The total torque includes two parts. One of the two
these approaches have also been developed to study parts, causing particle i to rotate, is contributed by the
granular flows in the transitional regime. An important tangential components of the traction distribution. Another,
method is the viscous, elastic-plastic continuum model. It often referred to as the rolling friction torque, is contributed
has been extensively applied to various cases (e.g., Wu, by the asymmetrical normal components of the traction
1990; Oda & Iwashita, 2000). distribution. It provides a resistance to relative rolling mo-
This paper analyses the three popular methods simu- tion between particles. Based on the forms of their stiffness
lating granular flows at a microscopic and macroscopic parameters, the force models can be divided into two
scales: discrete element method, averaging method and classes: linear model and nonlinear model. It is still open
viscous, elastic-plastic continuum model. The mathemati- for discussion which model is better. The detailed equa-
cal models of these methods and their applications to tions about the nonlinear model for spherical particles are
hopper flows are considered. The advantages and disad- listed in Table 1.
vantages involved in the theories and applications of these
Table 1 Equations for the calculation of interparticle force and torque
methods are discussed.
Force and torque Equations
2. Mathematical Models Normal elastic force fijne 4 *
Ei Rij* δ ( )
3
n 2
nij
ij
3
2.1 Discrete element method
)
1

In the DEM simulation, a granular material is modeled

Normal damping force fijnd
(
−cijn 8mij Ei* Rij*δ ijn 2
vijn

based on a finite number of discrete, semi-rigid spherical ⎛ 3

⎞ t
or polygon shaped particles interacting by means of con- Tangential elastic force fijte
⎜ ( (
− μ ij fijne ⎜ 1 − 1 − min vijt , δ ijmax ) δ ijmax ) 2
⎟ vˆ ij
⎟
⎝ ⎠
tact or non-contact forces, and the translational and rota-
( )
1
tional motions of every single particle in a considered sys- Tangential damping force fijtd 2c ijt 1.5 μij mij fijne 1 − vtij δ ijmax δ ijmax
2
vtij
tem are described by Newton’s laws of motion. For sim-
plicity, our present study is limited to coarse spherical par- Coulomb friction fijt − μ ij fijne vˆ ijt
ticle systems in which the effect of interstitial fluid and Friction torque m t
Cij × fijt
ij
non-contact forces such as the van der Waals and electro-
static forces can be ignored. Therefore, the governing Rolling friction torque mnij {
− min μ r ,ij fijne , μ r' ,ij ωijn ω
ˆ ijn }
equations for translational and rotational motion of particle i
where R = ( Ri R j ) ( Ri + R j ) , E = Ei ⎡⎣2(1 − ν )⎤⎦ , δ = Ri + R j − r j − ri ,
*
ij i
*
i
2 n
ij
can be given by
Vijn = ( Vij ⋅ nij ) ⋅ nij , Vijt = ( Vij × nij ) × nij , Vij = Vi − Vj + ωi × C ij − ω j × C ji ,
dv i
mi = ∑ fij + fib + mi g , (1)
dt j
nij = ( ri − r j ) ri − r j , vˆ tij = vtij vtij , ˆ nij = ωijn ωijn
ω ,

dω i δ ijmax = μ ij δ ijn (2 − ν i ) [2(1 − ν i )]

Ii = ∑ mij + mbi . (2)
dt j
To track the evolution of the positions of particles, their
The most complicated issue with this method is probably
initial positions and velocities are first given. The initial
the calculation of the interaction forces and torques be-
information is then used to calculate the interparticle forces
tween particles. Various approaches have been proposed
and torques for each contact by means of the force and
to model the interparticle forces. One of the most popular
torque models listed in Table 1. Solving Eqs. (1) and (2) by
force models is developed based on the consideration of
use of a numerical scheme gives new positions and ve-
contact elastic force and viscous contact damping force
(Langston et al., 1995a; Tanaka et al., 2002; Zhu & Yu, locities of these particles at time Δt. The time step, Δt,
2003). In general, the total interaction force between parti- should be set smaller than a certain critical value
cles can be expressed as a normal component and a tan- (Langston et al, 1995b). Using the new particle positions
gential one. The normal force includes elastic and damping and velocities, the calculation cycle is repeated in the next
components. The tangential force is described by static time step. Therefore, the evolution of the position of each
and dynamic frictions. Prior to the relative sliding of the particle can be determined.
contacting particles, the tangential force is described by the 2.2 Averaging method
static force including two components corresponding to the
normal elastic and damping forces respectively. When By use of a proper averaging procedure, the discrete
particles at contact start to slide relatively, the tangential system considered above can be transferred into a corre-
356 CHINA PARTICUOLOGY Vol. 3, No. 6, 2005

sponding continuum system. The balance equations of the macroscopic variables to the microscopic variables in the
continuum system have generally been derived by various discrete approach, and applied to the dynamical analysis
methods. Previous work often ignored the effect of the of granular flows (e.g., Langston et al., 1995a; Potapov &
rotational motion of particles. Therefore, the obtained bal- Campbell, 1996). However, some restrictive assumptions
ance equations only include those of mass and linear have to be employed in these models. Therefore, care
momentum. These equations are the same as those in the should be taken when using these models. For example, to
classical continuum mechanics, no matter which method a large degree, the volume averaging method proposed by
has been used to derive them. Recent studies have illus- Rothenburg and Selvadurai (1981) is only valid for
trated that in some cases, for example, shear band, the quasi-static systems because the inertial effect is ignored
gradient of particle rotation is very high (Oda & Iwashita, and a static equilibrium condition is used. Further, it is not
1999); consequently, some additional quantities should clear if the macroscopic properties obtained conform with
also be included in the continuum formulation to describe those in the balance equations in the continuum approach
this gradient. Based on such consideration, an extra equa- (i.e., Eqs. (3)–(5)). Similar problem also exists in the
tion has been derived to describe the rotation of the inde- time-volume averaging method of Zhang and Campbell
pendent particles. Therefore, complete balance equations (1992). The weighted time-volume averaging method
for a continuum system include those of mass, linear mo- proposed by Babic (1997) and further developed by Zhu
mentum and angular momentum. They correspond to the and Yu (2001 & 2002) (listed in Table 2) can overcome
balance equations of mass, translational motion and rota- these problems. According to the method, the DEM results
tional motion for individual particles, respectively. Accord- can be used to calculate the macroscopic variables such
ing to Babic (1997) and Zhu and Yu (2001 & 2002), these as mass density, velocity, stress and couple stress satis-
equations can be respectively given by fying the balance equations (3)–(5). However, this ap-
D( ρ ) + ρ∇ ⋅ u = 0 , (3) proach depends on the weighting function for small sys-
tems. The need to find out a suitable weighting function
D( ρ u) + ρ u∇ ⋅ u = ∇ ⋅ T + ρ g , (4)
has been noticed in our previous studies, and a weighting
D(λω) + λω∇ ⋅ u = ∇ ⋅ M + M' . (5) function has been recommended. The weighting function
The main macroscopic quantities involved in these bal- guarantees that the contribution of the particles nearer a
ance equations include mass density, velocity, angular probe point or time closer to a probe time to a considered
velocity, stress tensor and couple stress tensor. Various macroscopic property is larger (Zhu & Yu, 2002; 2003;
methods have been proposed to directly link these 2005b).

Table 2 Equations linking microscopic to macroscopic variables

Macroscopic variables* Equations linking microscopic to macroscopic variables

Mass density ρ ∫ ∑ hi mi ds
Tt i

1
Velocity u
ρ T∫
∑ h m v ds
i
i i i
t

1
Angular velocity ω
λ T∫
∑ h I ω ds
i
i i i
t

− ∫ ∑ hi mi v i '⊗ v i ' ds + ∫ ∑∑ g ij dij ⊗ fij ds + ∫ ∑ g ib dib ⊗ fib ds

Stress T Tt i Tt i j >i Tt i

1
− ∫ ∑ hi I i v i '⊗ ωi ds + ∑ ∑ gij dij ⊗ (mij − m ji )ds + ∫ ∑i gibdib ⊗ mib ds
Couple stress M
Tt i 2 T∫t i j >i Tt

1
Rate of supply of internal spin M ' ∑∑ (mij + m ji )(hi + hj )ds
2 T∫t i j >i
1 1
where hi = h(ri − r, s − t ) , λ = ∫ ∑ hi Ii ds , v i ' = v i − u , g ij = ∫ h( ri + rdij − r, s − t )dr , g ib = ∫ h(ri + rdib − r, s − t )dr
0 0
Tt i

*These macroscopic quantities depend on the weighting function h( r , t ) . A weighting function has been recom-
mended, and given by (Zhu and Yu, 2002; 2005b).
⎧ 1 c0Lt ⎡ 1⎛ L + t L + r ⎞⎤
⎪⎪ exp ⎢ − ⎜ ln2 t + ln2 p ⎟ ⎥ , ( r , t) ∈ Ω
h( r , t ) = ⎨ 2π 2π (L2t − t 2 )Lp (L2p − r 2 ) ⎢⎣ 2 ⎝⎜ Lt − t Lp − r ⎠⎟ ⎥⎦
⎪
⎪⎩ 0, others
 Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow 357

2.3 Viscous, elastic-plastic continuum model posed to derive the constitutive equations for granular
materials, i.e, the conventional plastic flow rule theory
Through the joint efforts of many scientists and engi-
(Rombach & Eibl, 1995; Wu & Schmidt, 1992) and the
neers, two types of constitutive models have been devel-
so-called double-shearing theory originated by Spencer
oped, namely plasticity flow models (e.g., Spencer, 1982;
(1964) and then extended by others (e.g., Hill & Wu, 1991,
Collins, 1990; Hill & Wu, 1992) and non-Newtonian fluid
1993b, 1996; Hill et al., 1999). The double shearing theory
models (Shen & Hopkins, 1988; Goldshtein & Shapiro,
assumes that the complete plastic deformation at any point
1995; Campbell, 1990). Moreover, a number of attempts
of the material consists of two shearing motions occurring
have been made to construct the so-called viscous, elas-
simultaneously on the α- and β-lines and which are su-
tic-plastic models by combining the plastic flow rule theory
with the elasticity theory and the fluid mechanics theory perposed one on the other, where the α- and β-lines are
(Rombach & Eibl, 1995; Schmidt & Wu, 1989). A common inclined at angles ±θα (θα = π 4 + φ 2 ) to the direction of
feature of this kind of models is that, the stress tensor is the principal stress σ1 in which φ is the angle of mate-
assumed to consist of a rate dependent part and a rate
independent part, namely rial internal friction. The plastic flow rule theory is based on
D D D a plastic potential gp, being a function of stresses, and
T = Tv + Ts , (6) assumes that, once the stress state in the material satisfies
D the yield condition f=0, plastic deformation occurs, and the
where T is the co-rotational rate of stress, given by
plastic strain rates d ijp are proportional to the derivatives
D ∂T
T= + u ⋅ ∇T + Tw − wT , (7) of the plastic potential with respect to the corresponding
∂t
stress, i.e.,
where w denotes the rate of spin tensor which is related to
the velocity vector u by the following geometric equations ∂g
dp = λ p . (13)
1 ∂T
w = (∇u − (∇u)T ) . (8)
2 The flow rule is called associated if gp is identical to the
D yield function f or non-associated otherwise. The most
The rate dependent part Tv is due to the viscous effect widely used yield condition is the so-called Mohr-Coulomb
and is determined either by the linear viscous fluid model failure criterion, which states that the material fails when
or a non-Newtonian fluid model. In the present work, we the shear stress τ acting on the surface of a material ele-
assume that the rate dependent part is linearly proportional ment attains the critical value
to the shear rate, similar to the deviatoric part of a Newto- τ = c + σ n tan φ , (14)
nian fluid, namely
D D
where σn is the compressive normal stress, and c is the
Tv = G d , (9) cohesion. Based on the Mohr-Coulomb criterion, various
where d is the rate of total deformation and the compo- constitutive equations have been developed utilising the
nents of G are given by associated flow rule (Cox et al., 1961). This kind of models
Gijrs = 2μ (δ irδ js − 31 δ ijrs ) , (10) usually predicts the stress field reasonably well. However,
unrealistically large dilation rates are usually predicted
in which μ is the viscosity of the material and δ denotes the
(Collins, 1990). Hence, two other types of models have
Dirac delta function.
been proposed. The first type of models abandon Cou-
The rate independent part of stress rate is assumed to
lomb's failure criterion and replace it by a family of yield
be generated by the material deformation rate d which
curves depending on the current density (Roscoe et al.,
consists of the elastic deformation rate component de and
1958). The second type of models use a plastic potential
the plastic deformation rate component dp. That is
different with the yield function (Lade, 1977; Schmidt & Wu,
d = de + dp . (11)
1989). Based on equations (11)–(14), the rate independent
The elastic deformation rate is related to the time rate of D
the stress tensor by Hooke’ law, so that part Ts can be determined by
dT D
de = D−1 , (12) Ts = Hd , (15)
dt
where D is the elasticity matrix, while the plastic deforma- where H is the so-called elastic-plastic matrix and its form
tion rate can be determined through the use of a frictional depends on the plasticity model used.
plasticity theory. The frictional plasticity theory is based on The other equations governing the velocity and stress
the assumption that the material flows according to a con- fields in granular media are the usual stress equations of
stitutive law if the stress state satisfies a yield condition. motion and the applied boundary conditions. In general,
Many different yield conditions have been established using the viscous, elastic-plastic continuum model, a typi-
(Farley & Valentin, 1968; Hill & Wu, 1993a; Matsuoka, cal granular flow problem can be described by the follow-
1984) and two different types of theories have been pro- ing boundary value problem,
358 CHINA PARTICUOLOGY Vol. 3, No. 6, 2005

BVP: find T and u such that microscopic or macroscopic scale (e.g., Wu, 1990;
Du Langston et al., 1995a; Xie & Shinohara, 1999; Zhu & Yu,
∇⋅T +F − ρ = 0 in Π, (16)
Dt 2004; 2005a; Steingart & Evans, 2005). Here, we just show
D D some of our results to illustrate the features of the ap-
T = Hd + G d in Π, (17) proaches introduced above.
u t =0 = u , T t =0 = T
o o
in Π, (18) 3.1 Discrete modelling
u = 0 on ∂Π a , T ⋅ n = f on ∂Π t , (19) We consider the steady state granular flow in a cylindri-
u ⋅ n = 0, (T ⋅ n) ⋅ t = − sgn(u ⋅ t )μ f ( T ⋅ n) ⋅ n on ∂Π f , (20) cal 3D hopper with flat bottom. Discrete simulation is per-
formed by means of DEM with the nonlinear force model.
where Π ∈ R , ∂Π = ∂Π a ∪ ∂Π t ∪ ∂Π f is the boundary of
2
The geometry of the hopper used in this work is shown in
Π with ∂Π a , ∂Π t and ∂Π f representing boundary with Fig. 1. It is cylindrical in shape, of diameter of 24d and with
a circular orifice of diameter of 8d at the centre of its flat
rigid constraints, prescribed traction and frictional force
bottom. 24,000 multi-sized spherical particles of uniform
respectively, n and t denote respectively the unit outward
size distribution in a range of 0.8 to 1.0d are used. In the
normal vector and unit tangential vector of ∂Π .
DEM simulation, particles with random initial velocities are
Due to the complexity of the underlying equations, so- first allowed to gradually settle onto the hopper under
phisticated numerical techniques are needed to solve gravity, giving a packing with a height of about 35d. These
boundary value problems of practical interest. Various particles are then discharged when the hopper outlet is
numerical schemes based on the Galerkin finite element opened. To achieve macroscopically steady state flow,
method have been developed to solve the static and dy-
particles discharged at the bottom outlet are recycled back
namic problems (Rombach & Eibl, 1995; Schmidt & Wu,
to the hopper randomly on the top of the particle bed.
1989). In comparison with the conventional finite element
These particles are set to the same velocities as the parti-
method for other mechanics disciplines, the new features
cles nearest to them in the hopper.
of the schemes for granular flows include the method
dealing with boundary friction, no-tension analysis and the 24 d
method dealing with both material and geometric nonlin-
earities.

3. Application to Hopper Flows 100 d

Granular flow in hoppers has been an important topic of

research worldwide for several decades (e.g., Polderman z
et al., 1987; Nedderman, 1995; Seville at el., 1997). The y
dynamic behavior of hopper flow is very complicated be- x
cause all the flow regimes about granular flows may coex-
ist. This flow has been studied extensively by means of 8d

both experimental technique and numerical simulation at a Fig. 1 Geometry of the hopper model used in the DEM simulation.

Scale: gd

(a) (b) (c)

Fig. 2 Internal flow and force structures when μp=0.6, cd=0.3, μw=0.3: (a) flow pattern; (b) velocity field; (c) normal force network. In Fig. 2(a),
particles charged in the hopper at different time slots are represented by different colours. In Fig. 2(b), each vector represents the com-
ponent of the velocity of a particle in a vertical mid-section as long as its centre is located between ± 1.0d of the section. In Fig. 2(c), the
thickness of each stick connecting two particle centres is proportional to the magnitude of the normal force between the two particles; the
thickness of the section is the same as that for Fig. 2(b).
 Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow 359

Since the trajectories of and contact forces acting on in- illustrative example, Fig. 3 shows the effects of the wall
dividual particles are traced in a DEM simulation, informa- friction coefficient, μ w , particle friction and damping coef-
tion for the flow pattern, velocity field and force structure of
ficients, μp and c d , on the statistical distribution of ver-
the hopper flow can be readily established, as shown in
Fig. 2. This indicates that DEM simulation can generate tical velocity. It can be observed that these mi-
microdynamic information which is difficult to obtain ex- cro-properties, especially, wall roughness, affect the flow
perimentally but key to elucidating the fundamentals gov- behaviour of particles. Detailed discusses can be seen in
erning the granular flow. For example, at present, the our previous work (Zhu & Yu, 2004).
forces among particles can only be measured under 3.2 Discrete-continuum transition
two-dimensional conditions or on the surface of a
three-dimensional flow (Løvoll et al., 1999; Geng et al., The values of the macroscopic quantities such as den-
2001). Such information can also be used to investigate sity, velocity, stress and couple stress can be obtained by
the physical parameters of hoppers and particles. As an use of the micro-dynamic information generated and the
above average method. Theoretically, the stress or couple
16
stress contains two parts: kinetic contribution related to the
transport of particles and collisional contribution related to
base values
12 μ
μ pp = 0.4
the interaction force or torque between particles and be-
tween particles and walls. For granular flow in a static re-
Probability density

c d = 0 .1
μμ ww = 0.1 gime, the kinetic contribution is very small and usually
8 ignored. However, for granular flow in a rapid flow regime,
the kinetic contribution is significant and can not be ignored.
To ensure the accuracy of the stress and couple stress
4
distribution, both of the two contributions are included in
the present calculation for the hopper flow in which the
three regimes may coexist. Since the resulting macro-
0
0 0.1 0.2 0.3 0.4 0.5 scopic variables should be reasonably axially symmetric,
Vertical velocity all computed macroscopic quantities are averaged from
Fig. 3 The dependence of vertical velocity on particle friction coeffi- the values of two symmetrical points and shown as a func-
cient, particle damping coefficient and wall friction coefficient tion of radius r and height z in cylindrical coordinates. The
(units for velocity are gd ). Base values are μp=0.6, cd=0.3, cylindrical coordinate {r , θ , z} starts at the center of the
and μw=0.3. orifice of the hopper.
Height z

Height z

Height z

Height z

Radius r Radius r Radius r Radius r

(a) (b) (c) (d)
Fig. 4 Contour plots of the distributions of: (a) vertical velocity uz; (b) mass density ρ ; (c) stress Tzz ; (d) couple stress M rθ when μp=0.6, cd=0.3, μw=0.3

(The units for velocity, mass density, stress, couple stress and radius/height are gd , πρp / 6 , πρpdg / 6 , πρp d 2 g / 6 and d, respectively).

The vertical velocity, uz , mass density, ρ , a repre- can be seen that these distributions are largely related to
the corresponding microscopic structures of the hopper
sentative stress, Tzz , and a representative couple stress,
flow. For example, since the particles near the orifice have
Mrθ , have been investigated in this cylindrical coordinates. large velocities, the averaged velocity is large in the orifice
The distributions of these properties are shown in Fig. 4. It region; corresponding to large normal forces occurring to
360 CHINA PARTICUOLOGY Vol. 3, No. 6, 2005

the particles around the bottom corner of the hopper, the maximum principal stress in the hopper area changes from
stress has a large magnitude in a region near the corner. a vertical direction to a direction nearly normal to the hop-
The stress or couple stress is not only related to the inter- per wall. Subsequently, the stresses decrease above the
action forces or torques, but also to the branch vectors outlet, and increase at the transition between the vertical
connecting the centers of mass of the two particles at wall and the inclined wall. Figure 7 shows the variation of
contact. Therefore, the macroscopic description is more wall pressure with time at some fixed spatial points during
complicated than the microscopic one. Similar to the mi- the process of material discharge. It is clear that the dis-
croscopic properties, the macroscopic properties should charge process starts with a wall pressure increase just
depend on the physical parameters of the hopper and above the outlet, then a decrease of pressure near the
particles. This is true as illustrated by Fig. 5, where the outlet occurs followed by a significant increase of pressure
effects of the wall roughness of the hopper, and the friction in the transition area. These results are consistent with
and damping coefficients between particles on the vertical reported experimental measurements (e.g., Blight, 1986).
normal stress are shown. Therefore, a continuum theory, to The results show that the model can capture the high
be general for granular flow, should be able to describe the pressure concentration in the transition area of a
effects of these parameters. bin-hopper. This provides important input for the optimal
design of hoppers to prevent possible structural failure
when in operation.

Fig. 5 The dependence of stress on particle friction coefficient, parti-

cle damping coefficient and wall friction coefficient at height
z = 9d (The units for the stresses and radius are respectively
πρpdg / 6 and d). Base values are μp=0.6, cd=0.3 and μw=0.3.

In the application of the combined approach of DEM and

averaging method, the main issue is that the number of
Fig. 6 Finite element mesh and principal stress distributions: (a) finite
particles which can be handled is limited by the present
element mesh; (b) stress at t=0.0 s; (c) stress at t=0.025 s.
computational capacity, and it is difficult to handle
non-spherical particles. Therefore, it is difficult to adapt this
Stations
approach to process modelling. Further, the obtained
macroscopic quantities depend on the selection of the size
of the cell where the discrete quantities are averaged. If the
size is too small, the resulting macroscopic quantities
fluctuate, giving unreasonable average results. On the
Pressure at stations

other hand, if it is too large, the local properties of these

quantities may be eliminated. The continuum approach
can overcome these problems, as illustrated below.
3.3 Continuum modelling
By use of the viscous, elastic-plastic continuum model
and a finite element method (FEM), we consider the flow of
granular solids through a silo with geometry as shown in
Fig. 6(a). Figures 6(b) and (c) show the principal stress
fields at two different instants. Among them, the stress field
Discharge time / μ s
at t=0 is the stress at the end of the filling condition, and
Fig. 7 Variation of wall pressure (kPa) with time (micro-second) at
also is the initial stress at the start of the discharge process.
several fixed spatial stations (see Fig. 6(a)).
During the process of discharge, the orientation of the
 Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow 361

p -1
4. Conclusions d plastic deformation rate component of d, m⋅s
-3
D elasticity matrix, N⋅m
Granular flow is a discrete system and, like a gas or liq- Ei -1 -2
modulus of elasticity of particle i, kg⋅m ⋅s
uid, can be described at different time and length scales. It
can be described by DEM on an individual particle or mi- fij interaction force between particles i and j, N
croscopic scale. By means of DEM simulation, detailed f i
b
interaction force between particle i and walls, N
information about the micro-dynamic behavior of the t
f tangential force between particle i and j, N
granular flow such as flow pattern, velocity, and flow and ij

ne
force structures can be obtained. These quantities are f ij normal elastic force between particle i and j, N
affected by the physical properties of particles and walls f yield function
such as friction and damping coefficients and wall rough- g gravitational acceleration, its magnitude, g, equal to
-2
ness. Granular material can also be described by contin- 9.81, m⋅s
uum model on a macroscopic scale. The macroscopic gp plastic potential
quantities such as velocity and stress can be obtained by G given by Eq. (10)
h weighting function
solving the balance equations by use of the computational -3
H elastic-plastic matrix, N⋅m
method such as FEM. However, this approach depends on 2
Ii moment of inertia, kg·m
the material constitutive equations used and cannot de-
scribe the effect of microscopic structure of granular flow. Lt , Lp parameters determining domain Ω
The combined approach of DEM simulation and averaging M couple stress tensor, Pa·m
method can overcome this problem. By use of the com- mi mass of particle i, kg
bined approach, the macro-dynamic behavior of granular 1
flow, including velocity, density, stress and couple stress mij =
2
( mi + m j )
can be determined. Similar to microscopic quantities, these
Mrθ component of couple stress tensor, Pa·m
macroscopic quantities obtained depend on the mi-
M' rate of supply of internal spin to particles, Pa
cro-properties of granular material. The combined ap-
mij torque acting on particle i exerted by particle j, N·m
proach takes into account the discrete nature of granular
materials and does not require any global assumption and mbi torque acting on particle i exerted by wall, N·m
thus allows a better understanding of the fundamental n unit outward normal vector of ∂Π
mechanisms of granular flow. The main disadvantage with r, ri position vector, m
this approach is that it is difficult to adapt it to process
r distance from probe point to the central axis of hopper, m
modelling because of the limited number of particles which
can be handled and the difficulty in handling non-spherical Ri radius of particles i, m
particles with current computer capacity. Further work is s time contributing to the averaged quantities, s
needed to develop an appropriate approach to overcome t time, s
these problems. Tt = [ −Lt + t , Lt + t ]
T stress tensor, Pa
Acknowledgement T
D
co-rotational rate of stress, Pa·s
-1

The authors are grateful to Australian Research Council and D D

Tv rate dependent part of T , Pa·s
-1
Chinese Academy of Sciences for the financial support of this work.
D D
Ts rate independent part of T , Pa·s
-1

Nomenclature Tzz vertical normal stress, Pa

c cohesion, Pa t unit tangential vector of ∂Π
c0 normalized constant of weighting function u velocity, m·s
-1

cd damping coefficient -1
uz component of mean velocity at vertical direction, m·s
cijn , cijt normal and tangential damping coefficients -1
vi velocity of particle i, m·s
Cij vector from the mass centre of the particle i to the
contact point with particle j, m v ijt vector of the relative tangential displacement of parti-
d maximum particle diameter, m cles i and j, m
-1
part of the branch vector connecting the mass centres w rate of spin tensor, 1·s
dij
z height of probe point, m
of particles i and j within Ω p , m
Greek letters
dbi ray from the mass center of particle i to ∂Ω p , via a ρ -3
mass density, kg·m
point on the wall and perpendicular to the tangential ρp particle density, kg·m
-3

plane of the point, m -1 -2

d rate of total deformation, m⋅s
-1 μ viscosity of material, kg·s ·m
d
e
elastic deformation rate component of d, m⋅s
-1 μw sliding friction coefficient between particles and walls
362 CHINA PARTICUOLOGY Vol. 3, No. 6, 2005

μp sliding friction coefficient between particles equations. J. Fluid Mech., 282, 75-114.
Herrmann, H. J. & Luding, S. (1998). Modeling granular media on
μ ij sliding friction coefficient between particles i and j the computer. Continuum Mech. Therm., 10, 189-231.
μr,ij Hill, J. M. & Wu, Y. H. (1991). Kinematically determined axially
rolling friction coefficient, m
symmetric plastic flows of metals and granular materials. Q. J.
μ '
r,ij rotational stiffness, m·N·s·rad
-1
Mech. Appl. Math., 44, 451-469.
-1 Hill, J. M. & Wu, Y. H. (1992). Some axially symmetric flows of
ωi angular velocity of particle i, rad·s
Mohr-Coulomb compressible granular materials. Proc. R. Soc.
ωnij component of the vector of the relative angular velocity Lond., A, 438, 67-93.
-1
of particles i and j in their contact plane, rad·s Hill, J. M. & Wu, Y. H. (1993a). Plastic flows of granular materials
Ω {
given by Ω = ( r , t ) r ∈ Ω p , t ∈ T } of shear index n: Part I yield functions. J. Mech. Phys. Solids, 41,
77-94.
Ωp domain in which quantities are averaged, given by Hill, J. M. & Wu, Y. H. (1993b). Plastic flows of granular materials
{
Ω p = r r ≤ Lp } of shear index n: Part II plane and axially symmetric problems
for n=2. J. Mech. Phys. Solids, 41, 95-115.
∂Ωp boundary of domain Ω p Hill, J. M. & Wu, Y. H. (1996). The punch problem for shear index
νi Poisson’s ratio of particle i granular materials. Q. J. Mech. Appl. Math., 49, 81-105.
Hill, J. M., Shi, J., Tordesillas, A. & Wu, Y. H. (1999). The velocity
σ1 principal stress, Pa field for the punch problem for dilatant granular materials. Q. J.
σn compressive normal stress, Pa Mech. Appl. Math., 52, 99-110.
τ shear stress, Pa Kishino, Y. (Ed.) (2001). Powders and Grains (a conference pro-
Π domain in R 2 ceeding containing numerous DEM studies on the physics of
granular matter). Lisse (The Netherland): A. A. Balkema Pub-
∂Π boundary of Π
lishers.
∂Π a boundary with rigid constraints Lade, P. V. (1977). Elasto-plastic stress strain theory for cohe-
∂Π t boundary with prescribed traction sionless soil with curved yield surface. Int. J. Solids Struct., 13,
1019-1035.
∂Π f boundary with frictional force Langston, P. A., Tüzün, U. & Heyes, D. M. (1995a). Discrete
φ angle of material internal friction element simulation of internal stress and flow fields in funnel
δ Dirac delta function flow hoppers. Powder Technol., 85, 153-169.
Langston, P. A., Tüzün, U. & Heyes, D. M. (1995b). Discrete
element simulation of granular flow in 2d and 3d hoppers: de-
References pendence of discharge rate and wall stress on particle interac-
Alonso-Marroquin, F. & Herrmann, H. J. (2005). Investigation of tions. Chem. Eng. Sci., 50, 967-987.
the incremental response of soils using a discrete element Løvoll, G., Måløy, K. J. & Flekkøy, E. G. (1999). Force measure-
model. J. Engng. Math., 52, 11-34. ments on static granular materials. Phys. Rev. E, 60,
Babic, M. (1997). Average balance equations for granular materi- 5872-5878.
als. Int. J. Eng. Sci., 35, 523-548. Lun, C. K. K., Savage, S. B., Jeffrey, D. J. & Chepurniy, N. (1984).
Blight, G. E. (1986). Pressures exerted by materials stored in silos: Kinetic theories for granular flow: inelastic particles in Couette
Part I coarse materials; Part II fine powders. Geotechnique, 36, flow and slightly inelastic particles in a general flowfield. J. Fluid
33-56. Mech., 140, 223-256.
Campbell, C. S. (1990). Rapid granular flows. Annu. Rev. Fluid Matsuoka, H. (1984). Deformation and strength of granular mate-
Mech., 22, 57-92. rials based on the theory of compounded mobilized planes and
Collins, J. (1990). Plane strain characteristics theory for soils and spatial mobilized plane. In Shahinpoor, M. (Ed.), Advances in
granular materials with density dependent yield criteria. J. Mech. the Mechanics and the Flow of Granular Materials (2, 814-836).
Phys. Solids, 38, 1-25. Houston (USA): Bulf Publishing Company.
Cox, A. D., Eason, G. & Hopkins, H. G. (1961). Axially symmetric Nedderman, R. M. (1992). Statics and Kinematics of Granular
plastic deformation in soils. Philos. Trans. R. Soc. Lond., A, 254, Materials. New York: Cambridge University Press.
1-45. Nedderman, R. M. (1995). The use of the kinematic model to
Cundall, P. A. & Strack, O. D. L. (1979). A discrete numerical predict the development of the stagnant zone boundary in the
model for granular assemblies. Geotechnique, 29, 47-65. batch discharge of a bunker. Chem. Eng. Sci., 50, 959-965.
de Gennes, P. G. (1999). Granular matter: a tentative view. Rev. Oda, M. & Iwashita, K. (1999). Mechanics of Granular Materials.
Mod. Phys., 71, S374-382. Rotterdam (The Netherland): A. A. Balkema Publishers.
Drescher, A. & de Josselin de Jong, G. (1972). Photoelastic veri- Oda, M. & Iwashita, K. (2000). Study on couple stress and shear
fication of a mechanical model for the flow of a granular material. band development in granular media based on numerical
J. Mech. Phys. Solids, 20, 337-351. simulation analyses. Int. J. Engng. Sci., 38, 1713-1740.
Farley, R. & Valentin, F. H. H. (1968). Effect of particle size upon Polderman, H. G., Boom, J., de Hilster, E. & Scott, A. M. (1987).
the strength of powders. Powder Technol., 1, 344-354. Solids flow velocity profiles in mass flow hoppers. Chem. Eng.
Geng, J. F., Longhi, E., Behringer, R. P. & Howell, D. W. (2001). Sci., 42, 737-744.
Memory in two-dimensional heap experiments. Phys. Rev. E, 64, Potapov, A. V. & Campbell, C. S. (1996). Computer simulation of
060301. hopper flow. Phys. Fluids, 8, 2884-2894.
Goldshtein, A. & Shapiro, M. (1995). Mechanics of collisional Ristow, G. H. (1992). Simulating granular flow with molecular
motion of granular materials, Part 1, general hydrodynamic dynamics. J. Phys. I, 2, 649-662.
 Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow 363

Rombach, G. & Eibl, J. (1995). Granular flow of materials in silos. Simulation of flow and mixing of particles in a rotating and
Int. J. Bulk Solid Handling, 15, 65-70. rocking cylinder. AIChE J., 44, 1266-1276.
Roscoe, K. H., Schofield, A. N. & Wroth, C. P. (1958). On the Wu, Y. H. (1990). Static and Dynamic Analyses of the Flow of Bulk
yielding of soils, Geotechnique, 9, 71-83. Materials Through Silos. PhD thesis, The University of Wollon-
Rothenburg, L. & Selvadurai, A. P. S. (1981). A micro-mechanical gong, Australia.
definition of the Cauchy stress tensor for particulate media. In Wu, Y. H. & Schmidt, L. C. (1992). A boundary element method for
Selvadurai, A. P. S. (Ed.), Proceedings of the International prediction of silo pressure. Int. J. Comput. Struct., 45, 315-323.
Symposium on the Mechanical Behavior of Structured Media Xie, H. Y. & Shinohara, K. (1999). Measurement of solids velocity
(pp.469-486). Ottawa, Canada. in a conical hopper by mass tracer particles. Chem. Eng. Sci.,
Schmidt, L. C. & Wu, Y. H. (1989). Prediction of dynamic wall 54, 455-459.
pressures on silos. Int. J. Bulk Solid Handling, 9, 333-338. Yamane, K., Nakagawa, M., Altobelli, S. A., Tanaka, T. & Tsuji, Y.
Seville, J. P. K., Tüzün, U. & Clift, R. (1997). Processing of Par- (1998). Steady particulate flows in a horizontal rotating cylinder.
ticulate Solids. Landon: Blackie Academic & Professional. Phys. Fluids, 10, 1419-1427.
Shen, H. & Hopkins, M. (1988). Stress in a rapid, simple shear Yang, R. Y., Zou, R. P. & Yu, A. B. (2003). Microdynamic analysis
flow of granular materials with multiple grain sizes. Part. Sci. of particle flow in a horizontal rotating drum. Powder Technol.,
Technol., 6, 1-15. 130, 138-146.
Spencer, A. J. M. (1964). A theory of the kinematics of ideal soils Zhang, Y. & Campbell, C. S. (1992). The interface between
under plane strain conditions. J. Mech. Phys. Solids, 12, fluid-like and solid-like behaviour in two-dimensional granular
337-351. flows. J. Fluid Mech., 237, 541-568.
Spencer, A. J. M. (1982). Deformation of ideal granular materials. Zhu, H. P. & Yu, A. B. (2001). A weighting function in the averag-
In Hopkins, H. G. & Sewell, M. J. (Eds.), Mechanics of Solids ing method of granular materials and its application to hopper
(pp. 607-652). New York: Pergamon Press. flow. Bulk Solids Handling, 21, 53-57.
Steingart, D. A. & Evans, J. W. (2005). Measurements of granular Zhu, H. P. & Yu, A. B. (2002). Averaging method of granular ma-
flows in two-dimensional hoppers by particle image velocimetry. terials. Phys. Rev. E, 66, 021302.
Part 1: experimental method and results. Chem. Eng. Sci., 60, Zhu, H. P. & Yu, A. B. (2003). The effects of wall and rolling re-
1043-1051. sistance on the couple stress of granular materials in vertical
Stewart, R. L., Bridgwater, J., Zhou, Y. C. & Yu, A. B. (2001). flow. Physica A, 325, 347-360.
Simulated and measured flow of granules in a bladed mixer ⎯ a Zhu, H. P. & Yu, A. B. (2004). Steady-state granular flow in a 3D
detailed comparison. Chem. Eng. Sci., 56, 5457-5471. cylindrical hopper with flat bottom using DEM simulation: mi-
Tanaka, K., Nishida, M., Kunimochi, T. & Takagi, T. (2002). Dis- croscopic analysis. J. Phys. D: Appl. Phys., 37, 1497-1508.
crete element simulation and experiment for dynamic response Zhu, H. P. & Yu, A. B. (2005a) Micromechanic analysis of hopper
of two-dismensional granular matter to the impact of a spherical flow of granular materials. J. Engng. Math., 52, 307-320.
projectile. Powder Technol., 124, 160-173. Zhu, H. P. & Yu, A. B. (2005b). Steady-state granular flow in a 3D
Walton, O. R. & Braun, R. L. (1986). Viscosity, granu- cylindrical hopper with flat bottom using DEM simulation: mac-
lar-temperature, and stress calculations for shearing assemblies roscopic analysis. Granular Matter, 7, 97-107.
of inelastic, frictional disks. J. Rheol., 30, 949-980.
Wightman, C., Moakher, M., Muzzio, F. J. & Walton, O. R. (1998). Manuscript received November 2, 2005 and accepted November 10, 2005.