- Email: [email protected]
Discrete and continuum modelling of granular flow
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: [email protected]
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 relationship 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 particles 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 The dynamic behavior of granular materials is very complicated. As with solids, they can withstand deformation and form heaps; as with liquids, they can flow; as with gases, they exhibit compressibility. These features give rise to another state of matter that is poorly understood (de Gennes, 1999). Corresponding to the fluid- and solid-like modes, different regimes have been identified in the past: quasi-static regime, rapid flow regime and a transitional regime that lies in between. To describe the flow behavior in these regimes, different approaches at different time and length scales have been developed in the past. Generally speaking, they can be categorize into two length scales: microscopic and macroscopic. At a microscopic or particle scale, granular material is a discrete system whose physical properties are discontinuous with respect to position and time. On the other hand, at a macroscopic or bulk scale, granular material is a continuum system whose physical properties are continuous. Based on such considerations, different approaches have been proposed to describe granular material at the two distinct scales. An important simulation approach which can describe granular flow at a particle scale is the so-called discrete element method (DEM) originally developed by Cundall and Strack (1979). The method considers a finite number of discrete interacting particles, and every particle in the considered system is described by Newton’s equations of motion related to translational and rotational motions. The motion of every particle is traced and the gradients of translation and rotation of a particle are determined in terms of the forces and the torques exerted on it.
DEM-based simulation has been used to study the physics of granular materials (e.g., Herrmann & Luding, 1998; Kishino, 2001; Zhu et al., 2004), and applied to investigate various engineering problems such as mixing (Wightman et al., 1998; Stewart et al., 2001), grinding/granulation (Yamane et al., 1998; Yang et al, 2003) and hopper flow (Ristow, 1992; Langston et al., 1995a; Zhu & Yu, 2004). By use of a proper averaging procedure, the discrete system considered above can be transferred into a corresponding continuum system. Extensive research has been carried out to develop such averaging methods. The methods thus far proposed can be divided into three classes: volume average (Drescher & de Josselin de Jong, 1972; Rothenburg & Selvadurai, 1981), time-volume average (Walton & Braun, 1986; Zhang & Campbell, 1992) and weighted time-volume average (Babic, 1997; Zhu & Yu, 2002). Based on a proper averaging method, the macroscopic quantities such as density, velocity and stress can be obtained in terms of the microscopic quantities such as velocities of particles, and interaction forces and torques between particles. The combined approach of DEM and averaging method has been applied to the macro-dynamical analysis of granular flows (e.g., Langston et al., 1995a; Potapov & Campbell, 1996; Zhu & Yu, 2005b), and has been used to study the intrinsic characteristics of granular materials such as the constitutive relationship under various flow conditions (Oda & Iwashita, 2000; Alonso-Marroquin & Herrmann, 2005). Granular material can also be modeled by continuum approach on a macroscopic scale. In the continuum approach, the macroscopic behavior of granular flow is described by the balance equations facilitated with constitu-
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow tive relations and boundary conditions. In the past, two continuum models developed based on the plasticity theory and kinetic theory of molecular dynamics have extensively been used to study the dynamic behavior of granular materials (Lun et al., 1984; Campbell, 1990; Nedderman, 1992). They have been illustrated to be applicable to quasi-static and rapid flow regimes, respectively. Recently, these approaches have also been developed to study granular flows in the transitional regime. An important method is the viscous, elastic-plastic continuum model. It has been extensively applied to various cases (e.g., Wu, 1990; Oda & Iwashita, 2000). This paper analyses the three popular methods simulating granular flows at a microscopic and macroscopic scales: discrete element method, averaging method and viscous, elastic-plastic continuum model. The mathematical models of these methods and their applications to hopper flows are considered. The advantages and disadvantages involved in the theories and applications of these methods are discussed.
355
force is described by the dynamic friction, usually, given by the Coulomb friction model. The contact between two particles is not a single point but is a finite area due to the deformation of both particles. The interparticle forces act over the contact region between particles rather than the mass centre of the particles, and they will generate a torque. The total torque includes two parts. One of the two parts, causing particle i to rotate, is contributed by the tangential components of the traction distribution. Another, often referred to as the rolling friction torque, is contributed by the asymmetrical normal components of the traction distribution. It provides a resistance to relative rolling motion between particles. Based on the forms of their stiffness parameters, the force models can be divided into two classes: linear model and nonlinear model. It is still open for discussion which model is better. The detailed equations about the nonlinear model for spherical particles are listed in Table 1. Table 1
Equations for the calculation of interparticle force and torque
Force and torque
2. Mathematical Models 2.1 Discrete element method In the DEM simulation, a granular material is modeled based on a finite number of discrete, semi-rigid spherical or polygon shaped particles interacting by means of contact or non-contact forces, and the translational and rotational motions of every single particle in a considered system are described by Newton’s laws of motion. For simplicity, our present study is limited to coarse spherical particle systems in which the effect of interstitial fluid and non-contact forces such as the van der Waals and electrostatic forces can be ignored. Therefore, the governing equations for translational and rotational motion of particle i can be given by dv i (1) mi = ∑ fij + fib + mi g , dt j Ii
dω i = ∑ mij + mbi . dt j
(2)
The most complicated issue with this method is probably the calculation of the interaction forces and torques between particles. Various approaches have been proposed to model the interparticle forces. One of the most popular force models is developed based on the consideration of contact elastic force and viscous contact damping force (Langston et al., 1995a; Tanaka et al., 2002; Zhu & Yu, 2003). In general, the total interaction force between particles can be expressed as a normal component and a tangential one. The normal force includes elastic and damping components. The tangential force is described by static and dynamic frictions. Prior to the relative sliding of the contacting particles, the tangential force is described by the static force including two components corresponding to the normal elastic and damping forces respectively. When particles at contact start to slide relatively, the tangential
Equations 3 n 2 ij
Normal elastic force fijne
4 * Ei Rij* δ 3
Normal damping force fijnd
−cijn 8mij Ei* Rij*δ ijn
Tangential elastic force fijte
⎛ − μ ij fijne ⎜ 1 − 1 − min vijt , δ ijmax ⎜ ⎝
Tangential damping force fijtd
2c ijt 1.5 μij mij fijne
Coulomb friction fijt
− μ ij fijne vˆ ijt
( )
nij
(
(
)
1 2
vijn
δ ijmax
)
3 2
1 − vtij δ ijmax δ ijmax
)
1 2
(
(
Friction torque m
Cij × fijt
Rolling friction torque mnij
ˆ ijn − min μ r ,ij fijne , μ r' ,ij ωijn ω
t ij
{
)
⎞ t ⎟ vˆ ij ⎟ ⎠ vtij
}
where R = ( Ri R j ) ( Ri + R j ) , E = Ei ⎡⎣2(1 − ν )⎤⎦ , δ = Ri + R j − r j − ri , Vijn = ( Vij ⋅ nij ) ⋅ nij , Vijt = ( Vij × nij ) × nij , Vij = Vi − Vj + ωi × C ij − ω j × C ji , * ij
nij = ( ri − r j ) ri − r j
* i
,
2 i
vˆ tij = vtij vtij
n ij
,
ˆ nij = ωijn ωijn ω
δ ijmax = μ ij δ ijn (2 − ν i ) [2(1 − ν i )]
To track the evolution of the positions of particles, their initial positions and velocities are first given. The initial information is then used to calculate the interparticle forces and torques for each contact by means of the force and torque models listed in Table 1. Solving Eqs. (1) and (2) by use of a numerical scheme gives new positions and velocities of these particles at time Δt. The time step, Δt, should be set smaller than a certain critical value (Langston et al, 1995b). Using the new particle positions and velocities, the calculation cycle is repeated in the next time step. Therefore, the evolution of the position of each particle can be determined.
2.2 Averaging method By use of a proper averaging procedure, the discrete system considered above can be transferred into a corre-
,
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
356
sponding continuum system. The balance equations of the continuum system have generally been derived by various methods. Previous work often ignored the effect of the rotational motion of particles. Therefore, the obtained balance equations only include those of mass and linear momentum. These equations are the same as those in the classical continuum mechanics, no matter which method has been used to derive them. Recent studies have illustrated that in some cases, for example, shear band, the gradient of particle rotation is very high (Oda & Iwashita, 1999); consequently, some additional quantities should also be included in the continuum formulation to describe this gradient. Based on such consideration, an extra equation has been derived to describe the rotation of the independent particles. Therefore, complete balance equations for a continuum system include those of mass, linear momentum and angular momentum. They correspond to the balance equations of mass, translational motion and rotational motion for individual particles, respectively. According to Babic (1997) and Zhu and Yu (2001 & 2002), these equations can be respectively given by (3) D( ρ ) + ρ∇ ⋅ u = 0 ,
D( ρ u) + ρ u∇ ⋅ u = ∇ ⋅ T + ρ g ,
(4)
D(λω) + λω∇ ⋅ u = ∇ ⋅ M + M' .
(5)
The main macroscopic quantities involved in these balance equations include mass density, velocity, angular velocity, stress tensor and couple stress tensor. Various methods have been proposed to directly link these Table 2
macroscopic variables to the microscopic variables in the discrete approach, and applied to the dynamical analysis of granular flows (e.g., Langston et al., 1995a; Potapov & Campbell, 1996). However, some restrictive assumptions have to be employed in these models. Therefore, care should be taken when using these models. For example, to a large degree, the volume averaging method proposed by Rothenburg and Selvadurai (1981) is only valid for quasi-static systems because the inertial effect is ignored and a static equilibrium condition is used. Further, it is not clear if the macroscopic properties obtained conform with those in the balance equations in the continuum approach (i.e., Eqs. (3)–(5)). Similar problem also exists in the time-volume averaging method of Zhang and Campbell (1992). The weighted time-volume averaging method proposed by Babic (1997) and further developed by Zhu and Yu (2001 & 2002) (listed in Table 2) can overcome these problems. According to the method, the DEM results can be used to calculate the macroscopic variables such as mass density, velocity, stress and couple stress satisfying the balance equations (3)–(5). However, this approach depends on the weighting function for small systems. The need to find out a suitable weighting function has been noticed in our previous studies, and a weighting function has been recommended. The weighting function guarantees that the contribution of the particles nearer a probe point or time closer to a probe time to a considered macroscopic property is larger (Zhu & Yu, 2002; 2003; 2005b).
Equations linking microscopic to macroscopic variables
Macroscopic variables*
Equations linking microscopic to macroscopic variables
∫ ∑ hi mi ds
Mass density ρ
i
Tt
1
Velocity u
ρ T∫
∑ h m v ds
1
∑ h I ω ds
i
λ T∫
Angular velocity ω
i
i
i
t
i i
i
i
t
Stress T
− ∫ ∑ hi mi v i '⊗ v i ' ds + ∫ ∑∑ g ij dij ⊗ fij ds + ∫ ∑ g ib dib ⊗ fib ds
Couple stress M
− ∫ ∑ hi I i v i '⊗ ωi ds +
i
Tt
Tt
Rate of supply of internal spin M '
i
Tt
j >i
i
Tt
i
1 ∑ ∑ gij dij ⊗ (mij − m ji )ds + ∫ ∑i gibdib ⊗ mib ds 2 T∫t i j >i Tt
1 ∑∑ (mij + m ji )(hi + hj )ds 2 T∫t i j >i 1
1
0
0
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 Tt
i
*These macroscopic quantities depend on the weighting function h( r , t ) . A weighting function has been recommended, and given by (Zhu and Yu, 2002; 2005b). ⎧ 1 ⎡ 1⎛ L + t L + r ⎞⎤ c0Lt exp ⎢ − ⎜ ln2 t + ln2 p ⎟ ⎥ , ( r , t) ∈ Ω ⎪⎪ h( r , t ) = ⎨ 2π 2π (L2t − t 2 )Lp (L2p − r 2 ) Lp − r ⎠⎟ ⎥⎦ ⎢⎣ 2 ⎝⎜ Lt − t ⎪ others ⎪⎩ 0,
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow
2.3 Viscous, elastic-plastic continuum model Through the joint efforts of many scientists and engineers, two types of constitutive models have been developed, namely plasticity flow models (e.g., Spencer, 1982; Collins, 1990; Hill & Wu, 1992) and non-Newtonian fluid models (Shen & Hopkins, 1988; Goldshtein & Shapiro, 1995; Campbell, 1990). Moreover, a number of attempts have been made to construct the so-called viscous, elastic-plastic models by combining the plastic flow rule theory with the elasticity theory and the fluid mechanics theory (Rombach & Eibl, 1995; Schmidt & Wu, 1989). A common feature of this kind of models is that, the stress tensor is assumed to consist of a rate dependent part and a rate independent part, namely D
D
D
T = Tv + Ts ,
(6)
D
where T is the co-rotational rate of stress, given by D ∂T (7) T= + u ⋅ ∇T + Tw − wT , ∂t where w denotes the rate of spin tensor which is related to the velocity vector u by the following geometric equations 1 (8) w = (∇u − (∇u)T ) . 2 D
The rate dependent part Tv is due to the viscous effect and is determined either by the linear viscous fluid model or a non-Newtonian fluid model. In the present work, we assume that the rate dependent part is linearly proportional to the shear rate, similar to the deviatoric part of a Newtonian fluid, namely D
D
Tv = G d , (9) where d is the rate of total deformation and the components of G are given by (10) Gijrs = 2μ (δ irδ js − 31 δ ijrs ) , in which μ is the viscosity of the material and δ denotes the Dirac delta function. The rate independent part of stress rate is assumed to be generated by the material deformation rate d which consists of the elastic deformation rate component de and the plastic deformation rate component dp. That is d = de + dp . (11) The elastic deformation rate is related to the time rate of the stress tensor by Hooke’ law, so that dT , (12) de = D−1 dt where D is the elasticity matrix, while the plastic deformation rate can be determined through the use of a frictional plasticity theory. The frictional plasticity theory is based on the assumption that the material flows according to a constitutive law if the stress state satisfies a yield condition. Many different yield conditions have been established (Farley & Valentin, 1968; Hill & Wu, 1993a; Matsuoka, 1984) and two different types of theories have been pro-
357
posed to derive the constitutive equations for granular materials, i.e, the conventional plastic flow rule theory (Rombach & Eibl, 1995; Wu & Schmidt, 1992) and the so-called double-shearing theory originated by Spencer (1964) and then extended by others (e.g., Hill & Wu, 1991, 1993b, 1996; Hill et al., 1999). The double shearing theory assumes that the complete plastic deformation at any point of the material consists of two shearing motions occurring simultaneously on the α- and β-lines and which are superposed one on the other, where the α- and β-lines are inclined at angles ±θα (θα = π 4 + φ 2 ) to the direction of the principal stress σ1 in which φ is the angle of material internal friction. The plastic flow rule theory is based on a plastic potential gp, being a function of stresses, and assumes that, once the stress state in the material satisfies the yield condition f=0, plastic deformation occurs, and the plastic strain rates d ijp are proportional to the derivatives of the plastic potential with respect to the corresponding stress, i.e., ∂g dp = λ p . (13) ∂T The flow rule is called associated if gp is identical to the yield function f or non-associated otherwise. The most widely used yield condition is the so-called Mohr-Coulomb failure criterion, which states that the material fails when the shear stress τ acting on the surface of a material element attains the critical value τ = c + σ n tan φ , (14) where σn is the compressive normal stress, and c is the cohesion. Based on the Mohr-Coulomb criterion, various constitutive equations have been developed utilising the associated flow rule (Cox et al., 1961). This kind of models usually predicts the stress field reasonably well. However, unrealistically large dilation rates are usually predicted (Collins, 1990). Hence, two other types of models have been proposed. The first type of models abandon Coulomb's failure criterion and replace it by a family of yield curves depending on the current density (Roscoe et al., 1958). The second type of models use a plastic potential different with the yield function (Lade, 1977; Schmidt & Wu, 1989). Based on equations (11)–(14), the rate independent D
part Ts can be determined by D
Ts = Hd ,
(15)
where H is the so-called elastic-plastic matrix and its form depends on the plasticity model used. The other equations governing the velocity and stress fields in granular media are the usual stress equations of motion and the applied boundary conditions. In general, using the viscous, elastic-plastic continuum model, a typical granular flow problem can be described by the following boundary value problem,
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
358
BVP: find T and u such that Du ∇⋅T +F − ρ = 0 in Π, Dt D
D
u t =0 = u , T t =0 = T
o
(17) (18)
3.1 Discrete modelling
(19)
We consider the steady state granular flow in a cylindrical 3D hopper with flat bottom. Discrete simulation is performed by means of DEM with the nonlinear force model. The geometry of the hopper used in this work is shown in 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 bottom. 24,000 multi-sized spherical particles of uniform size distribution in a range of 0.8 to 1.0d are used. In the DEM simulation, particles with random initial velocities are first allowed to gradually settle onto the hopper under gravity, giving a packing with a height of about 35d. These particles are then discharged when the hopper outlet is opened. To achieve macroscopically steady state flow, particles discharged at the bottom outlet are recycled back to the hopper randomly on the top of the particle bed. These particles are set to the same velocities as the particles nearest to them in the hopper.
(16)
T = Hd + G d in Π, o
microscopic or macroscopic scale (e.g., Wu, 1990; Langston et al., 1995a; Xie & Shinohara, 1999; Zhu & Yu, 2004; 2005a; Steingart & Evans, 2005). Here, we just show some of our results to illustrate the features of the approaches introduced above.
in Π,
u = 0 on ∂Π a , T ⋅ n = f on ∂Π t , u ⋅ n = 0, (T ⋅ n) ⋅ t = − sgn(u ⋅ t )μ f ( T ⋅ n) ⋅ n on ∂Π f ,
(20)
where Π ∈ R , ∂Π = ∂Π a ∪ ∂Π t ∪ ∂Π f is the boundary of 2
Π with ∂Π a , ∂Π t and ∂Π f representing boundary with rigid constraints, prescribed traction and frictional force respectively, n and t denote respectively the unit outward normal vector and unit tangential vector of ∂Π . Due to the complexity of the underlying equations, sophisticated numerical techniques are needed to solve boundary value problems of practical interest. Various numerical schemes based on the Galerkin finite element method have been developed to solve the static and dynamic problems (Rombach & Eibl, 1995; Schmidt & Wu, 1989). In comparison with the conventional finite element method for other mechanics disciplines, the new features of the schemes for granular flows include the method dealing with boundary friction, no-tension analysis and the method dealing with both material and geometric nonlinearities.
24 d
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 et al., 1987; Nedderman, 1995; Seville at el., 1997). The dynamic behavior of hopper flow is very complicated because all the flow regimes about granular flows may coexist. This flow has been studied extensively by means of both experimental technique and numerical simulation at a Scale:
(a)
Fig. 2
z y x 8d
Fig. 1
Geometry of the hopper model used in the DEM simulation.
gd
(b)
(c)
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 component 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 Since the trajectories of and contact forces acting on individual particles are traced in a DEM simulation, information for the flow pattern, velocity field and force structure of the hopper flow can be readily established, as shown in Fig. 2. This indicates that DEM simulation can generate microdynamic information which is difficult to obtain experimentally but key to elucidating the fundamentals governing the granular flow. For example, at present, the forces among particles can only be measured under two-dimensional conditions or on the surface of a three-dimensional flow (Løvoll et al., 1999; Geng et al., 2001). Such information can also be used to investigate the physical parameters of hoppers and particles. As an 16
base values
μ μ pp = 0.4
c d = 0 .1
μμ ww = 0.1 8
4
0 0
0.1
0.2
0.4
0.3
0.5
Vertical velocity
Fig. 3
The dependence of vertical velocity on particle friction coefficient, particle damping coefficient and wall friction coefficient (units for velocity are
gd ). Base values are μp=0.6, cd=0.3,
and μw=0.3.
ficients, μp and c d , on the statistical distribution of vertical velocity. It can be observed that these micro-properties, especially, wall roughness, affect the flow behaviour of particles. Detailed discusses can be seen in our previous work (Zhu & Yu, 2004).
3.2 Discrete-continuum transition The values of the macroscopic quantities such as density, velocity, stress and couple stress can be obtained by use of the micro-dynamic information generated and the above average method. Theoretically, the stress or couple stress contains two parts: kinetic contribution related to the transport of particles and collisional contribution related to the interaction force or torque between particles and between particles and walls. For granular flow in a static regime, the kinetic contribution is very small and usually 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 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 macroscopic variables should be reasonably axially symmetric, all computed macroscopic quantities are averaged from the values of two symmetrical points and shown as a function of radius r and height z in cylindrical coordinates. The cylindrical coordinate {r , θ , z} starts at the center of the
Height z
Height z
orifice of the hopper.
Height z Fig. 4
illustrative example, Fig. 3 shows the effects of the wall friction coefficient, μ w , particle friction and damping coef-
Height z
Probability density
12
359
Radius r
Radius r
Radius r
Radius r
(a)
(b)
(c)
(d)
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
The vertical velocity, uz , mass density, ρ , a representative stress, Tzz , and a representative couple stress, Mrθ , have been investigated in this cylindrical coordinates.
The distributions of these properties are shown in Fig. 4. It
gd , πρp / 6 , πρpdg / 6 , πρp d 2 g / 6 and d, respectively).
can be seen that these distributions are largely related to the corresponding microscopic structures of the hopper flow. For example, since the particles near the orifice have large velocities, the averaged velocity is large in the orifice region; corresponding to large normal forces occurring to
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
360
the particles around the bottom corner of the hopper, the stress has a large magnitude in a region near the corner. The stress or couple stress is not only related to the interaction forces or torques, but also to the branch vectors connecting the centers of mass of the two particles at contact. Therefore, the macroscopic description is more complicated than the microscopic one. Similar to the microscopic properties, the macroscopic properties should depend on the physical parameters of the hopper and particles. This is true as illustrated by Fig. 5, where the effects of the wall roughness of the hopper, and the friction and damping coefficients between particles on the vertical normal stress are shown. Therefore, a continuum theory, to be general for granular flow, should be able to describe the effects of these parameters.
Fig. 5
maximum principal stress in the hopper area changes from a vertical direction to a direction nearly normal to the hopper wall. Subsequently, the stresses decrease above the outlet, and increase at the transition between the vertical wall and the inclined wall. Figure 7 shows the variation of wall pressure with time at some fixed spatial points during the process of material discharge. It is clear that the discharge process starts with a wall pressure increase just above the outlet, then a decrease of pressure near the outlet occurs followed by a significant increase of pressure in the transition area. These results are consistent with reported experimental measurements (e.g., Blight, 1986). The results show that the model can capture the high pressure concentration in the transition area of a bin-hopper. This provides important input for the optimal design of hoppers to prevent possible structural failure when in operation.
The dependence of stress on particle friction coefficient, particle 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.
Fig. 6
Stations
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 at t=0 is the stress at the end of the filling condition, and also is the initial stress at the start of the discharge process. During the process of discharge, the orientation of the
Finite element mesh and principal stress distributions: (a) finite element mesh; (b) stress at t=0.0 s; (c) stress at t=0.025 s.
Pressure at stations
In the application of the combined approach of DEM and averaging method, the main issue is that the number of particles which can be handled is limited by the present computational capacity, and it is difficult to handle non-spherical particles. Therefore, it is difficult to adapt this 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 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.
Discharge time / μ s
Fig. 7
Variation of wall pressure (kPa) with time (micro-second) at several fixed spatial stations (see Fig. 6(a)).
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow
4. Conclusions Granular flow is a discrete system and, like a gas or liquid, can be described at different time and length scales. It can be described by DEM on an individual particle or microscopic scale. By means of DEM simulation, detailed information about the micro-dynamic behavior of the granular flow such as flow pattern, velocity, and flow and force structures can be obtained. These quantities are affected by the physical properties of particles and walls such as friction and damping coefficients and wall roughness. Granular material can also be described by continuum model on a macroscopic scale. The macroscopic quantities such as velocity and stress can be obtained by solving the balance equations by use of the computational method such as FEM. However, this approach depends on the material constitutive equations used and cannot describe the effect of microscopic structure of granular flow. The combined approach of DEM simulation and averaging method can overcome this problem. By use of the combined approach, the macro-dynamic behavior of granular flow, including velocity, density, stress and couple stress can be determined. Similar to microscopic quantities, these macroscopic quantities obtained depend on the micro-properties of granular material. The combined 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. The main disadvantage with this approach is that it is difficult to adapt it to process modelling because of the limited number of particles which can be handled and the difficulty in handling non-spherical particles with current computer capacity. Further work is needed to develop an appropriate approach to overcome these problems.
Acknowledgement The authors are grateful to Australian Research Council and Chinese Academy of Sciences for the financial support of this work.
p
361
d D
plastic deformation rate component of d, m⋅s -3 elasticity matrix, N⋅m
Ei
modulus of elasticity of particle i, kg⋅m ⋅s
fij
interaction force between particles i and j, N
-1
c c0 cd
cohesion, Pa normalized constant of weighting function damping coefficient
cijn , cijt
normal and tangential damping coefficients
Cij
vector from the mass centre of the particle i to the contact point with particle j, m maximum particle diameter, m part of the branch vector connecting the mass centres of particles i and j within Ω p , m
d dij dbi
ray from the mass center of particle i to ∂Ω p , via a
d e d
point on the wall and perpendicular to the tangential plane of the point, m -1 rate of total deformation, m⋅s -1 elastic deformation rate component of d, m⋅s
-2
b i
interaction force between particle i and walls, N
t ij
tangential force between particle i and j, N
ne ij
f
normal elastic force between particle i and j, N
f g gp G h H
yield function gravitational acceleration, its magnitude, g, equal to -2 9.81, m⋅s plastic potential given by Eq. (10) weighting function -3 elastic-plastic matrix, N⋅m
Ii
moment of inertia, kg·m
f
f
2
Lt , Lp
parameters determining domain Ω
M mi
couple stress tensor, Pa·m
mij
=
mass of particle i, kg 1 ( mi + m j ) 2
Mrθ
component of couple stress tensor, Pa·m
M' mij
rate of supply of internal spin to particles, Pa
mbi
torque acting on particle i exerted by wall, N·m
torque acting on particle i exerted by particle j, N·m
n
unit outward normal vector of ∂Π
r, ri
position vector, m
r
distance from probe point to the central axis of hopper, m
Ri
radius of particles i, m
s
time contributing to the averaged quantities, s
t
time, s
Tt
= [ −Lt + t , Lt + t ]
T
stress tensor, Pa
D
T
co-rotational rate of stress, Pa·s
D
Tv
D
rate dependent part of T , Pa·s
D
Nomenclature
-1
D
-1
-1
Ts
rate independent part of T , Pa·s
Tzz
vertical normal stress, Pa
t u
unit tangential vector of ∂Π -1 velocity, m·s
uz
component of mean velocity at vertical direction, m·s
vi
velocity of particle i, m·s
v ijt
vector of the relative tangential displacement of particles i and j, m -1 rate of spin tensor, 1·s height of probe point, m
w z
-1
-1
-1
Greek letters ρ
mass density, kg·m
ρp
particle density, kg·m
μ μw
viscosity of material, kg·s ·m
-3 -3 -1
-2
sliding friction coefficient between particles and walls
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
362 μp
sliding friction coefficient between particles
μ ij
sliding friction coefficient between particles i and j
μr,ij
rolling friction coefficient, m
μ
rotational stiffness, m·N·s·rad
' r,ij
-1 -1
ωi
angular velocity of particle i, rad·s
ωnij
component of the vector of the relative angular velocity -1 of particles i and j in their contact plane, rad·s
Ω
given by Ω = ( r , t ) r ∈ Ω p , t ∈ T
Ωp
domain in which quantities are averaged, given by
{
{
Ω p = r r ≤ Lp
}
}
∂Ωp
boundary of domain Ω p
νi
Poisson’s ratio of particle i
σ1
principal stress, Pa
σn τ Π
compressive normal stress, Pa shear stress, Pa
∂Π
domain in R 2 boundary of Π
∂Π a
boundary with rigid constraints
∂Π t
boundary with prescribed traction
∂Π f
boundary with frictional force
φ
angle of material internal friction Dirac delta function
δ
References Alonso-Marroquin, F. & Herrmann, H. J. (2005). Investigation of the incremental response of soils using a discrete element model. J. Engng. Math., 52, 11-34. Babic, M. (1997). Average balance equations for granular materials. Int. J. Eng. Sci., 35, 523-548. Blight, G. E. (1986). Pressures exerted by materials stored in silos: Part I coarse materials; Part II fine powders. Geotechnique, 36, 33-56. Campbell, C. S. (1990). Rapid granular flows. Annu. Rev. Fluid Mech., 22, 57-92. Collins, J. (1990). Plane strain characteristics theory for soils and granular materials with density dependent yield criteria. J. Mech. Phys. Solids, 38, 1-25. Cox, A. D., Eason, G. & Hopkins, H. G. (1961). Axially symmetric plastic deformation in soils. Philos. Trans. R. Soc. Lond., A, 254, 1-45. Cundall, P. A. & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29, 47-65. de Gennes, P. G. (1999). Granular matter: a tentative view. Rev. Mod. Phys., 71, S374-382. Drescher, A. & de Josselin de Jong, G. (1972). Photoelastic verification of a mechanical model for the flow of a granular material. J. Mech. Phys. Solids, 20, 337-351. Farley, R. & Valentin, F. H. H. (1968). Effect of particle size upon the strength of powders. Powder Technol., 1, 344-354. Geng, J. F., Longhi, E., Behringer, R. P. & Howell, D. W. (2001). Memory in two-dimensional heap experiments. Phys. Rev. E, 64, 060301. Goldshtein, A. & Shapiro, M. (1995). Mechanics of collisional motion of granular materials, Part 1, general hydrodynamic
equations. J. Fluid Mech., 282, 75-114. Herrmann, H. J. & Luding, S. (1998). Modeling granular media on the computer. Continuum Mech. Therm., 10, 189-231. Hill, J. M. & Wu, Y. H. (1991). Kinematically determined axially symmetric plastic flows of metals and granular materials. Q. J. Mech. Appl. Math., 44, 451-469. Hill, J. M. & Wu, Y. H. (1992). Some axially symmetric flows of Mohr-Coulomb compressible granular materials. Proc. R. Soc. Lond., A, 438, 67-93. Hill, J. M. & Wu, Y. H. (1993a). Plastic flows of granular materials of shear index n: Part I yield functions. J. Mech. Phys. Solids, 41, 77-94. Hill, J. M. & Wu, Y. H. (1993b). Plastic flows of granular materials of shear index n: Part II plane and axially symmetric problems for n=2. J. Mech. Phys. Solids, 41, 95-115. Hill, J. M. & Wu, Y. H. (1996). The punch problem for shear index granular materials. Q. J. Mech. Appl. Math., 49, 81-105. Hill, J. M., Shi, J., Tordesillas, A. & Wu, Y. H. (1999). The velocity field for the punch problem for dilatant granular materials. Q. J. Mech. Appl. Math., 52, 99-110. Kishino, Y. (Ed.) (2001). Powders and Grains (a conference proceeding containing numerous DEM studies on the physics of granular matter). Lisse (The Netherland): A. A. Balkema Publishers. Lade, P. V. (1977). Elasto-plastic stress strain theory for cohesionless soil with curved yield surface. Int. J. Solids Struct., 13, 1019-1035. Langston, P. A., Tüzün, U. & Heyes, D. M. (1995a). Discrete element simulation of internal stress and flow fields in funnel 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: dependence of discharge rate and wall stress on particle interactions. Chem. Eng. Sci., 50, 967-987. Løvoll, G., Måløy, K. J. & Flekkøy, E. G. (1999). Force measurements on static granular materials. Phys. Rev. E, 60, 5872-5878. Lun, C. K. K., Savage, S. B., Jeffrey, D. J. & Chepurniy, N. (1984). Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J. Fluid Mech., 140, 223-256. Matsuoka, H. (1984). Deformation and strength of granular materials based on the theory of compounded mobilized planes and spatial mobilized plane. In Shahinpoor, M. (Ed.), Advances in the Mechanics and the Flow of Granular Materials (2, 814-836). Houston (USA): Bulf Publishing Company. Nedderman, R. M. (1992). Statics and Kinematics of Granular Materials. New York: Cambridge University Press. Nedderman, R. M. (1995). The use of the kinematic model to predict the development of the stagnant zone boundary in the batch discharge of a bunker. Chem. Eng. Sci., 50, 959-965. Oda, M. & Iwashita, K. (1999). Mechanics of Granular Materials. Rotterdam (The Netherland): A. A. Balkema Publishers. Oda, M. & Iwashita, K. (2000). Study on couple stress and shear band development in granular media based on numerical simulation analyses. Int. J. Engng. Sci., 38, 1713-1740. Polderman, H. G., Boom, J., de Hilster, E. & Scott, A. M. (1987). Solids flow velocity profiles in mass flow hoppers. Chem. Eng. Sci., 42, 737-744. Potapov, A. V. & Campbell, C. S. (1996). Computer simulation of hopper flow. Phys. Fluids, 8, 2884-2894. Ristow, G. H. (1992). Simulating granular flow with molecular dynamics. J. Phys. I, 2, 649-662.
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow Rombach, G. & Eibl, J. (1995). Granular flow of materials in silos. Int. J. Bulk Solid Handling, 15, 65-70. Roscoe, K. H., Schofield, A. N. & Wroth, C. P. (1958). On the yielding of soils, Geotechnique, 9, 71-83. Rothenburg, L. & Selvadurai, A. P. S. (1981). A micro-mechanical definition of the Cauchy stress tensor for particulate media. In Selvadurai, A. P. S. (Ed.), Proceedings of the International Symposium on the Mechanical Behavior of Structured Media (pp.469-486). Ottawa, Canada. Schmidt, L. C. & Wu, Y. H. (1989). Prediction of dynamic wall pressures on silos. Int. J. Bulk Solid Handling, 9, 333-338. Seville, J. P. K., Tüzün, U. & Clift, R. (1997). Processing of Particulate Solids. Landon: Blackie Academic & Professional. Shen, H. & Hopkins, M. (1988). Stress in a rapid, simple shear flow of granular materials with multiple grain sizes. Part. Sci. Technol., 6, 1-15. Spencer, A. J. M. (1964). A theory of the kinematics of ideal soils under plane strain conditions. J. Mech. Phys. Solids, 12, 337-351. Spencer, A. J. M. (1982). Deformation of ideal granular materials. In Hopkins, H. G. & Sewell, M. J. (Eds.), Mechanics of Solids (pp. 607-652). New York: Pergamon Press. Steingart, D. A. & Evans, J. W. (2005). Measurements of granular flows in two-dimensional hoppers by particle image velocimetry. Part 1: experimental method and results. Chem. Eng. Sci., 60, 1043-1051. Stewart, R. L., Bridgwater, J., Zhou, Y. C. & Yu, A. B. (2001). Simulated and measured flow of granules in a bladed mixer ⎯ a detailed comparison. Chem. Eng. Sci., 56, 5457-5471. Tanaka, K., Nishida, M., Kunimochi, T. & Takagi, T. (2002). Discrete element simulation and experiment for dynamic response of two-dismensional granular matter to the impact of a spherical projectile. Powder Technol., 124, 160-173. Walton, O. R. & Braun, R. L. (1986). Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheol., 30, 949-980. Wightman, C., Moakher, M., Muzzio, F. J. & Walton, O. R. (1998).
363
Simulation of flow and mixing of particles in a rotating and rocking cylinder. AIChE J., 44, 1266-1276. Wu, Y. H. (1990). Static and Dynamic Analyses of the Flow of Bulk Materials Through Silos. PhD thesis, The University of Wollongong, Australia. Wu, Y. H. & Schmidt, L. C. (1992). A boundary element method for prediction of silo pressure. Int. J. Comput. Struct., 45, 315-323. Xie, H. Y. & Shinohara, K. (1999). Measurement of solids velocity in a conical hopper by mass tracer particles. Chem. Eng. Sci., 54, 455-459. Yamane, K., Nakagawa, M., Altobelli, S. A., Tanaka, T. & Tsuji, Y. (1998). Steady particulate flows in a horizontal rotating cylinder. Phys. Fluids, 10, 1419-1427. Yang, R. Y., Zou, R. P. & Yu, A. B. (2003). Microdynamic analysis of particle flow in a horizontal rotating drum. Powder Technol., 130, 138-146. Zhang, Y. & Campbell, C. S. (1992). The interface between fluid-like and solid-like behaviour in two-dimensional granular flows. J. Fluid Mech., 237, 541-568. Zhu, H. P. & Yu, A. B. (2001). A weighting function in the averaging method of granular materials and its application to hopper flow. Bulk Solids Handling, 21, 53-57. Zhu, H. P. & Yu, A. B. (2002). Averaging method of granular materials. Phys. Rev. E, 66, 021302. Zhu, H. P. & Yu, A. B. (2003). The effects of wall and rolling resistance on the couple stress of granular materials in vertical flow. Physica A, 325, 347-360. Zhu, H. P. & Yu, A. B. (2004). Steady-state granular flow in a 3D cylindrical hopper with flat bottom using DEM simulation: microscopic analysis. J. Phys. D: Appl. Phys., 37, 1497-1508. Zhu, H. P. & Yu, A. B. (2005a) Micromechanic analysis of hopper flow of granular materials. J. Engng. Math., 52, 307-320. Zhu, H. P. & Yu, A. B. (2005b). Steady-state granular flow in a 3D cylindrical hopper with flat bottom using DEM simulation: macroscopic analysis. Granular Matter, 7, 97-107. Manuscript received November 2, 2005 and accepted November 10, 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: [email protected]
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 relationship 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 particles 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 The dynamic behavior of granular materials is very complicated. As with solids, they can withstand deformation and form heaps; as with liquids, they can flow; as with gases, they exhibit compressibility. These features give rise to another state of matter that is poorly understood (de Gennes, 1999). Corresponding to the fluid- and solid-like modes, different regimes have been identified in the past: quasi-static regime, rapid flow regime and a transitional regime that lies in between. To describe the flow behavior in these regimes, different approaches at different time and length scales have been developed in the past. Generally speaking, they can be categorize into two length scales: microscopic and macroscopic. At a microscopic or particle scale, granular material is a discrete system whose physical properties are discontinuous with respect to position and time. On the other hand, at a macroscopic or bulk scale, granular material is a continuum system whose physical properties are continuous. Based on such considerations, different approaches have been proposed to describe granular material at the two distinct scales. An important simulation approach which can describe granular flow at a particle scale is the so-called discrete element method (DEM) originally developed by Cundall and Strack (1979). The method considers a finite number of discrete interacting particles, and every particle in the considered system is described by Newton’s equations of motion related to translational and rotational motions. The motion of every particle is traced and the gradients of translation and rotation of a particle are determined in terms of the forces and the torques exerted on it.
DEM-based simulation has been used to study the physics of granular materials (e.g., Herrmann & Luding, 1998; Kishino, 2001; Zhu et al., 2004), and applied to investigate various engineering problems such as mixing (Wightman et al., 1998; Stewart et al., 2001), grinding/granulation (Yamane et al., 1998; Yang et al, 2003) and hopper flow (Ristow, 1992; Langston et al., 1995a; Zhu & Yu, 2004). By use of a proper averaging procedure, the discrete system considered above can be transferred into a corresponding continuum system. Extensive research has been carried out to develop such averaging methods. The methods thus far proposed can be divided into three classes: volume average (Drescher & de Josselin de Jong, 1972; Rothenburg & Selvadurai, 1981), time-volume average (Walton & Braun, 1986; Zhang & Campbell, 1992) and weighted time-volume average (Babic, 1997; Zhu & Yu, 2002). Based on a proper averaging method, the macroscopic quantities such as density, velocity and stress can be obtained in terms of the microscopic quantities such as velocities of particles, and interaction forces and torques between particles. The combined approach of DEM and averaging method has been applied to the macro-dynamical analysis of granular flows (e.g., Langston et al., 1995a; Potapov & Campbell, 1996; Zhu & Yu, 2005b), and has been used to study the intrinsic characteristics of granular materials such as the constitutive relationship under various flow conditions (Oda & Iwashita, 2000; Alonso-Marroquin & Herrmann, 2005). Granular material can also be modeled by continuum approach on a macroscopic scale. In the continuum approach, the macroscopic behavior of granular flow is described by the balance equations facilitated with constitu-
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow tive relations and boundary conditions. In the past, two continuum models developed based on the plasticity theory and kinetic theory of molecular dynamics have extensively been used to study the dynamic behavior of granular materials (Lun et al., 1984; Campbell, 1990; Nedderman, 1992). They have been illustrated to be applicable to quasi-static and rapid flow regimes, respectively. Recently, these approaches have also been developed to study granular flows in the transitional regime. An important method is the viscous, elastic-plastic continuum model. It has been extensively applied to various cases (e.g., Wu, 1990; Oda & Iwashita, 2000). This paper analyses the three popular methods simulating granular flows at a microscopic and macroscopic scales: discrete element method, averaging method and viscous, elastic-plastic continuum model. The mathematical models of these methods and their applications to hopper flows are considered. The advantages and disadvantages involved in the theories and applications of these methods are discussed.
355
force is described by the dynamic friction, usually, given by the Coulomb friction model. The contact between two particles is not a single point but is a finite area due to the deformation of both particles. The interparticle forces act over the contact region between particles rather than the mass centre of the particles, and they will generate a torque. The total torque includes two parts. One of the two parts, causing particle i to rotate, is contributed by the tangential components of the traction distribution. Another, often referred to as the rolling friction torque, is contributed by the asymmetrical normal components of the traction distribution. It provides a resistance to relative rolling motion between particles. Based on the forms of their stiffness parameters, the force models can be divided into two classes: linear model and nonlinear model. It is still open for discussion which model is better. The detailed equations about the nonlinear model for spherical particles are listed in Table 1. Table 1
Equations for the calculation of interparticle force and torque
Force and torque
2. Mathematical Models 2.1 Discrete element method In the DEM simulation, a granular material is modeled based on a finite number of discrete, semi-rigid spherical or polygon shaped particles interacting by means of contact or non-contact forces, and the translational and rotational motions of every single particle in a considered system are described by Newton’s laws of motion. For simplicity, our present study is limited to coarse spherical particle systems in which the effect of interstitial fluid and non-contact forces such as the van der Waals and electrostatic forces can be ignored. Therefore, the governing equations for translational and rotational motion of particle i can be given by dv i (1) mi = ∑ fij + fib + mi g , dt j Ii
dω i = ∑ mij + mbi . dt j
(2)
The most complicated issue with this method is probably the calculation of the interaction forces and torques between particles. Various approaches have been proposed to model the interparticle forces. One of the most popular force models is developed based on the consideration of contact elastic force and viscous contact damping force (Langston et al., 1995a; Tanaka et al., 2002; Zhu & Yu, 2003). In general, the total interaction force between particles can be expressed as a normal component and a tangential one. The normal force includes elastic and damping components. The tangential force is described by static and dynamic frictions. Prior to the relative sliding of the contacting particles, the tangential force is described by the static force including two components corresponding to the normal elastic and damping forces respectively. When particles at contact start to slide relatively, the tangential
Equations 3 n 2 ij
Normal elastic force fijne
4 * Ei Rij* δ 3
Normal damping force fijnd
−cijn 8mij Ei* Rij*δ ijn
Tangential elastic force fijte
⎛ − μ ij fijne ⎜ 1 − 1 − min vijt , δ ijmax ⎜ ⎝
Tangential damping force fijtd
2c ijt 1.5 μij mij fijne
Coulomb friction fijt
− μ ij fijne vˆ ijt
( )
nij
(
(
)
1 2
vijn
δ ijmax
)
3 2
1 − vtij δ ijmax δ ijmax
)
1 2
(
(
Friction torque m
Cij × fijt
Rolling friction torque mnij
ˆ ijn − min μ r ,ij fijne , μ r' ,ij ωijn ω
t ij
{
)
⎞ t ⎟ vˆ ij ⎟ ⎠ vtij
}
where R = ( Ri R j ) ( Ri + R j ) , E = Ei ⎡⎣2(1 − ν )⎤⎦ , δ = Ri + R j − r j − ri , Vijn = ( Vij ⋅ nij ) ⋅ nij , Vijt = ( Vij × nij ) × nij , Vij = Vi − Vj + ωi × C ij − ω j × C ji , * ij
nij = ( ri − r j ) ri − r j
* i
,
2 i
vˆ tij = vtij vtij
n ij
,
ˆ nij = ωijn ωijn ω
δ ijmax = μ ij δ ijn (2 − ν i ) [2(1 − ν i )]
To track the evolution of the positions of particles, their initial positions and velocities are first given. The initial information is then used to calculate the interparticle forces and torques for each contact by means of the force and torque models listed in Table 1. Solving Eqs. (1) and (2) by use of a numerical scheme gives new positions and velocities of these particles at time Δt. The time step, Δt, should be set smaller than a certain critical value (Langston et al, 1995b). Using the new particle positions and velocities, the calculation cycle is repeated in the next time step. Therefore, the evolution of the position of each particle can be determined.
2.2 Averaging method By use of a proper averaging procedure, the discrete system considered above can be transferred into a corre-
,
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
356
sponding continuum system. The balance equations of the continuum system have generally been derived by various methods. Previous work often ignored the effect of the rotational motion of particles. Therefore, the obtained balance equations only include those of mass and linear momentum. These equations are the same as those in the classical continuum mechanics, no matter which method has been used to derive them. Recent studies have illustrated that in some cases, for example, shear band, the gradient of particle rotation is very high (Oda & Iwashita, 1999); consequently, some additional quantities should also be included in the continuum formulation to describe this gradient. Based on such consideration, an extra equation has been derived to describe the rotation of the independent particles. Therefore, complete balance equations for a continuum system include those of mass, linear momentum and angular momentum. They correspond to the balance equations of mass, translational motion and rotational motion for individual particles, respectively. According to Babic (1997) and Zhu and Yu (2001 & 2002), these equations can be respectively given by (3) D( ρ ) + ρ∇ ⋅ u = 0 ,
D( ρ u) + ρ u∇ ⋅ u = ∇ ⋅ T + ρ g ,
(4)
D(λω) + λω∇ ⋅ u = ∇ ⋅ M + M' .
(5)
The main macroscopic quantities involved in these balance equations include mass density, velocity, angular velocity, stress tensor and couple stress tensor. Various methods have been proposed to directly link these Table 2
macroscopic variables to the microscopic variables in the discrete approach, and applied to the dynamical analysis of granular flows (e.g., Langston et al., 1995a; Potapov & Campbell, 1996). However, some restrictive assumptions have to be employed in these models. Therefore, care should be taken when using these models. For example, to a large degree, the volume averaging method proposed by Rothenburg and Selvadurai (1981) is only valid for quasi-static systems because the inertial effect is ignored and a static equilibrium condition is used. Further, it is not clear if the macroscopic properties obtained conform with those in the balance equations in the continuum approach (i.e., Eqs. (3)–(5)). Similar problem also exists in the time-volume averaging method of Zhang and Campbell (1992). The weighted time-volume averaging method proposed by Babic (1997) and further developed by Zhu and Yu (2001 & 2002) (listed in Table 2) can overcome these problems. According to the method, the DEM results can be used to calculate the macroscopic variables such as mass density, velocity, stress and couple stress satisfying the balance equations (3)–(5). However, this approach depends on the weighting function for small systems. The need to find out a suitable weighting function has been noticed in our previous studies, and a weighting function has been recommended. The weighting function guarantees that the contribution of the particles nearer a probe point or time closer to a probe time to a considered macroscopic property is larger (Zhu & Yu, 2002; 2003; 2005b).
Equations linking microscopic to macroscopic variables
Macroscopic variables*
Equations linking microscopic to macroscopic variables
∫ ∑ hi mi ds
Mass density ρ
i
Tt
1
Velocity u
ρ T∫
∑ h m v ds
1
∑ h I ω ds
i
λ T∫
Angular velocity ω
i
i
i
t
i i
i
i
t
Stress T
− ∫ ∑ hi mi v i '⊗ v i ' ds + ∫ ∑∑ g ij dij ⊗ fij ds + ∫ ∑ g ib dib ⊗ fib ds
Couple stress M
− ∫ ∑ hi I i v i '⊗ ωi ds +
i
Tt
Tt
Rate of supply of internal spin M '
i
Tt
j >i
i
Tt
i
1 ∑ ∑ gij dij ⊗ (mij − m ji )ds + ∫ ∑i gibdib ⊗ mib ds 2 T∫t i j >i Tt
1 ∑∑ (mij + m ji )(hi + hj )ds 2 T∫t i j >i 1
1
0
0
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 Tt
i
*These macroscopic quantities depend on the weighting function h( r , t ) . A weighting function has been recommended, and given by (Zhu and Yu, 2002; 2005b). ⎧ 1 ⎡ 1⎛ L + t L + r ⎞⎤ c0Lt exp ⎢ − ⎜ ln2 t + ln2 p ⎟ ⎥ , ( r , t) ∈ Ω ⎪⎪ h( r , t ) = ⎨ 2π 2π (L2t − t 2 )Lp (L2p − r 2 ) Lp − r ⎠⎟ ⎥⎦ ⎢⎣ 2 ⎝⎜ Lt − t ⎪ others ⎪⎩ 0,
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow
2.3 Viscous, elastic-plastic continuum model Through the joint efforts of many scientists and engineers, two types of constitutive models have been developed, namely plasticity flow models (e.g., Spencer, 1982; Collins, 1990; Hill & Wu, 1992) and non-Newtonian fluid models (Shen & Hopkins, 1988; Goldshtein & Shapiro, 1995; Campbell, 1990). Moreover, a number of attempts have been made to construct the so-called viscous, elastic-plastic models by combining the plastic flow rule theory with the elasticity theory and the fluid mechanics theory (Rombach & Eibl, 1995; Schmidt & Wu, 1989). A common feature of this kind of models is that, the stress tensor is assumed to consist of a rate dependent part and a rate independent part, namely D
D
D
T = Tv + Ts ,
(6)
D
where T is the co-rotational rate of stress, given by D ∂T (7) T= + u ⋅ ∇T + Tw − wT , ∂t where w denotes the rate of spin tensor which is related to the velocity vector u by the following geometric equations 1 (8) w = (∇u − (∇u)T ) . 2 D
The rate dependent part Tv is due to the viscous effect and is determined either by the linear viscous fluid model or a non-Newtonian fluid model. In the present work, we assume that the rate dependent part is linearly proportional to the shear rate, similar to the deviatoric part of a Newtonian fluid, namely D
D
Tv = G d , (9) where d is the rate of total deformation and the components of G are given by (10) Gijrs = 2μ (δ irδ js − 31 δ ijrs ) , in which μ is the viscosity of the material and δ denotes the Dirac delta function. The rate independent part of stress rate is assumed to be generated by the material deformation rate d which consists of the elastic deformation rate component de and the plastic deformation rate component dp. That is d = de + dp . (11) The elastic deformation rate is related to the time rate of the stress tensor by Hooke’ law, so that dT , (12) de = D−1 dt where D is the elasticity matrix, while the plastic deformation rate can be determined through the use of a frictional plasticity theory. The frictional plasticity theory is based on the assumption that the material flows according to a constitutive law if the stress state satisfies a yield condition. Many different yield conditions have been established (Farley & Valentin, 1968; Hill & Wu, 1993a; Matsuoka, 1984) and two different types of theories have been pro-
357
posed to derive the constitutive equations for granular materials, i.e, the conventional plastic flow rule theory (Rombach & Eibl, 1995; Wu & Schmidt, 1992) and the so-called double-shearing theory originated by Spencer (1964) and then extended by others (e.g., Hill & Wu, 1991, 1993b, 1996; Hill et al., 1999). The double shearing theory assumes that the complete plastic deformation at any point of the material consists of two shearing motions occurring simultaneously on the α- and β-lines and which are superposed one on the other, where the α- and β-lines are inclined at angles ±θα (θα = π 4 + φ 2 ) to the direction of the principal stress σ1 in which φ is the angle of material internal friction. The plastic flow rule theory is based on a plastic potential gp, being a function of stresses, and assumes that, once the stress state in the material satisfies the yield condition f=0, plastic deformation occurs, and the plastic strain rates d ijp are proportional to the derivatives of the plastic potential with respect to the corresponding stress, i.e., ∂g dp = λ p . (13) ∂T The flow rule is called associated if gp is identical to the yield function f or non-associated otherwise. The most widely used yield condition is the so-called Mohr-Coulomb failure criterion, which states that the material fails when the shear stress τ acting on the surface of a material element attains the critical value τ = c + σ n tan φ , (14) where σn is the compressive normal stress, and c is the cohesion. Based on the Mohr-Coulomb criterion, various constitutive equations have been developed utilising the associated flow rule (Cox et al., 1961). This kind of models usually predicts the stress field reasonably well. However, unrealistically large dilation rates are usually predicted (Collins, 1990). Hence, two other types of models have been proposed. The first type of models abandon Coulomb's failure criterion and replace it by a family of yield curves depending on the current density (Roscoe et al., 1958). The second type of models use a plastic potential different with the yield function (Lade, 1977; Schmidt & Wu, 1989). Based on equations (11)–(14), the rate independent D
part Ts can be determined by D
Ts = Hd ,
(15)
where H is the so-called elastic-plastic matrix and its form depends on the plasticity model used. The other equations governing the velocity and stress fields in granular media are the usual stress equations of motion and the applied boundary conditions. In general, using the viscous, elastic-plastic continuum model, a typical granular flow problem can be described by the following boundary value problem,
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
358
BVP: find T and u such that Du ∇⋅T +F − ρ = 0 in Π, Dt D
D
u t =0 = u , T t =0 = T
o
(17) (18)
3.1 Discrete modelling
(19)
We consider the steady state granular flow in a cylindrical 3D hopper with flat bottom. Discrete simulation is performed by means of DEM with the nonlinear force model. The geometry of the hopper used in this work is shown in 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 bottom. 24,000 multi-sized spherical particles of uniform size distribution in a range of 0.8 to 1.0d are used. In the DEM simulation, particles with random initial velocities are first allowed to gradually settle onto the hopper under gravity, giving a packing with a height of about 35d. These particles are then discharged when the hopper outlet is opened. To achieve macroscopically steady state flow, particles discharged at the bottom outlet are recycled back to the hopper randomly on the top of the particle bed. These particles are set to the same velocities as the particles nearest to them in the hopper.
(16)
T = Hd + G d in Π, o
microscopic or macroscopic scale (e.g., Wu, 1990; Langston et al., 1995a; Xie & Shinohara, 1999; Zhu & Yu, 2004; 2005a; Steingart & Evans, 2005). Here, we just show some of our results to illustrate the features of the approaches introduced above.
in Π,
u = 0 on ∂Π a , T ⋅ n = f on ∂Π t , u ⋅ n = 0, (T ⋅ n) ⋅ t = − sgn(u ⋅ t )μ f ( T ⋅ n) ⋅ n on ∂Π f ,
(20)
where Π ∈ R , ∂Π = ∂Π a ∪ ∂Π t ∪ ∂Π f is the boundary of 2
Π with ∂Π a , ∂Π t and ∂Π f representing boundary with rigid constraints, prescribed traction and frictional force respectively, n and t denote respectively the unit outward normal vector and unit tangential vector of ∂Π . Due to the complexity of the underlying equations, sophisticated numerical techniques are needed to solve boundary value problems of practical interest. Various numerical schemes based on the Galerkin finite element method have been developed to solve the static and dynamic problems (Rombach & Eibl, 1995; Schmidt & Wu, 1989). In comparison with the conventional finite element method for other mechanics disciplines, the new features of the schemes for granular flows include the method dealing with boundary friction, no-tension analysis and the method dealing with both material and geometric nonlinearities.
24 d
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 et al., 1987; Nedderman, 1995; Seville at el., 1997). The dynamic behavior of hopper flow is very complicated because all the flow regimes about granular flows may coexist. This flow has been studied extensively by means of both experimental technique and numerical simulation at a Scale:
(a)
Fig. 2
z y x 8d
Fig. 1
Geometry of the hopper model used in the DEM simulation.
gd
(b)
(c)
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 component 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 Since the trajectories of and contact forces acting on individual particles are traced in a DEM simulation, information for the flow pattern, velocity field and force structure of the hopper flow can be readily established, as shown in Fig. 2. This indicates that DEM simulation can generate microdynamic information which is difficult to obtain experimentally but key to elucidating the fundamentals governing the granular flow. For example, at present, the forces among particles can only be measured under two-dimensional conditions or on the surface of a three-dimensional flow (Løvoll et al., 1999; Geng et al., 2001). Such information can also be used to investigate the physical parameters of hoppers and particles. As an 16
base values
μ μ pp = 0.4
c d = 0 .1
μμ ww = 0.1 8
4
0 0
0.1
0.2
0.4
0.3
0.5
Vertical velocity
Fig. 3
The dependence of vertical velocity on particle friction coefficient, particle damping coefficient and wall friction coefficient (units for velocity are
gd ). Base values are μp=0.6, cd=0.3,
and μw=0.3.
ficients, μp and c d , on the statistical distribution of vertical velocity. It can be observed that these micro-properties, especially, wall roughness, affect the flow behaviour of particles. Detailed discusses can be seen in our previous work (Zhu & Yu, 2004).
3.2 Discrete-continuum transition The values of the macroscopic quantities such as density, velocity, stress and couple stress can be obtained by use of the micro-dynamic information generated and the above average method. Theoretically, the stress or couple stress contains two parts: kinetic contribution related to the transport of particles and collisional contribution related to the interaction force or torque between particles and between particles and walls. For granular flow in a static regime, the kinetic contribution is very small and usually 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 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 macroscopic variables should be reasonably axially symmetric, all computed macroscopic quantities are averaged from the values of two symmetrical points and shown as a function of radius r and height z in cylindrical coordinates. The cylindrical coordinate {r , θ , z} starts at the center of the
Height z
Height z
orifice of the hopper.
Height z Fig. 4
illustrative example, Fig. 3 shows the effects of the wall friction coefficient, μ w , particle friction and damping coef-
Height z
Probability density
12
359
Radius r
Radius r
Radius r
Radius r
(a)
(b)
(c)
(d)
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
The vertical velocity, uz , mass density, ρ , a representative stress, Tzz , and a representative couple stress, Mrθ , have been investigated in this cylindrical coordinates.
The distributions of these properties are shown in Fig. 4. It
gd , πρp / 6 , πρpdg / 6 , πρp d 2 g / 6 and d, respectively).
can be seen that these distributions are largely related to the corresponding microscopic structures of the hopper flow. For example, since the particles near the orifice have large velocities, the averaged velocity is large in the orifice region; corresponding to large normal forces occurring to
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
360
the particles around the bottom corner of the hopper, the stress has a large magnitude in a region near the corner. The stress or couple stress is not only related to the interaction forces or torques, but also to the branch vectors connecting the centers of mass of the two particles at contact. Therefore, the macroscopic description is more complicated than the microscopic one. Similar to the microscopic properties, the macroscopic properties should depend on the physical parameters of the hopper and particles. This is true as illustrated by Fig. 5, where the effects of the wall roughness of the hopper, and the friction and damping coefficients between particles on the vertical normal stress are shown. Therefore, a continuum theory, to be general for granular flow, should be able to describe the effects of these parameters.
Fig. 5
maximum principal stress in the hopper area changes from a vertical direction to a direction nearly normal to the hopper wall. Subsequently, the stresses decrease above the outlet, and increase at the transition between the vertical wall and the inclined wall. Figure 7 shows the variation of wall pressure with time at some fixed spatial points during the process of material discharge. It is clear that the discharge process starts with a wall pressure increase just above the outlet, then a decrease of pressure near the outlet occurs followed by a significant increase of pressure in the transition area. These results are consistent with reported experimental measurements (e.g., Blight, 1986). The results show that the model can capture the high pressure concentration in the transition area of a bin-hopper. This provides important input for the optimal design of hoppers to prevent possible structural failure when in operation.
The dependence of stress on particle friction coefficient, particle 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.
Fig. 6
Stations
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 at t=0 is the stress at the end of the filling condition, and also is the initial stress at the start of the discharge process. During the process of discharge, the orientation of the
Finite element mesh and principal stress distributions: (a) finite element mesh; (b) stress at t=0.0 s; (c) stress at t=0.025 s.
Pressure at stations
In the application of the combined approach of DEM and averaging method, the main issue is that the number of particles which can be handled is limited by the present computational capacity, and it is difficult to handle non-spherical particles. Therefore, it is difficult to adapt this 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 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.
Discharge time / μ s
Fig. 7
Variation of wall pressure (kPa) with time (micro-second) at several fixed spatial stations (see Fig. 6(a)).
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow
4. Conclusions Granular flow is a discrete system and, like a gas or liquid, can be described at different time and length scales. It can be described by DEM on an individual particle or microscopic scale. By means of DEM simulation, detailed information about the micro-dynamic behavior of the granular flow such as flow pattern, velocity, and flow and force structures can be obtained. These quantities are affected by the physical properties of particles and walls such as friction and damping coefficients and wall roughness. Granular material can also be described by continuum model on a macroscopic scale. The macroscopic quantities such as velocity and stress can be obtained by solving the balance equations by use of the computational method such as FEM. However, this approach depends on the material constitutive equations used and cannot describe the effect of microscopic structure of granular flow. The combined approach of DEM simulation and averaging method can overcome this problem. By use of the combined approach, the macro-dynamic behavior of granular flow, including velocity, density, stress and couple stress can be determined. Similar to microscopic quantities, these macroscopic quantities obtained depend on the micro-properties of granular material. The combined 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. The main disadvantage with this approach is that it is difficult to adapt it to process modelling because of the limited number of particles which can be handled and the difficulty in handling non-spherical particles with current computer capacity. Further work is needed to develop an appropriate approach to overcome these problems.
Acknowledgement The authors are grateful to Australian Research Council and Chinese Academy of Sciences for the financial support of this work.
p
361
d D
plastic deformation rate component of d, m⋅s -3 elasticity matrix, N⋅m
Ei
modulus of elasticity of particle i, kg⋅m ⋅s
fij
interaction force between particles i and j, N
-1
c c0 cd
cohesion, Pa normalized constant of weighting function damping coefficient
cijn , cijt
normal and tangential damping coefficients
Cij
vector from the mass centre of the particle i to the contact point with particle j, m maximum particle diameter, m part of the branch vector connecting the mass centres of particles i and j within Ω p , m
d dij dbi
ray from the mass center of particle i to ∂Ω p , via a
d e d
point on the wall and perpendicular to the tangential plane of the point, m -1 rate of total deformation, m⋅s -1 elastic deformation rate component of d, m⋅s
-2
b i
interaction force between particle i and walls, N
t ij
tangential force between particle i and j, N
ne ij
f
normal elastic force between particle i and j, N
f g gp G h H
yield function gravitational acceleration, its magnitude, g, equal to -2 9.81, m⋅s plastic potential given by Eq. (10) weighting function -3 elastic-plastic matrix, N⋅m
Ii
moment of inertia, kg·m
f
f
2
Lt , Lp
parameters determining domain Ω
M mi
couple stress tensor, Pa·m
mij
=
mass of particle i, kg 1 ( mi + m j ) 2
Mrθ
component of couple stress tensor, Pa·m
M' mij
rate of supply of internal spin to particles, Pa
mbi
torque acting on particle i exerted by wall, N·m
torque acting on particle i exerted by particle j, N·m
n
unit outward normal vector of ∂Π
r, ri
position vector, m
r
distance from probe point to the central axis of hopper, m
Ri
radius of particles i, m
s
time contributing to the averaged quantities, s
t
time, s
Tt
= [ −Lt + t , Lt + t ]
T
stress tensor, Pa
D
T
co-rotational rate of stress, Pa·s
D
Tv
D
rate dependent part of T , Pa·s
D
Nomenclature
-1
D
-1
-1
Ts
rate independent part of T , Pa·s
Tzz
vertical normal stress, Pa
t u
unit tangential vector of ∂Π -1 velocity, m·s
uz
component of mean velocity at vertical direction, m·s
vi
velocity of particle i, m·s
v ijt
vector of the relative tangential displacement of particles i and j, m -1 rate of spin tensor, 1·s height of probe point, m
w z
-1
-1
-1
Greek letters ρ
mass density, kg·m
ρp
particle density, kg·m
μ μw
viscosity of material, kg·s ·m
-3 -3 -1
-2
sliding friction coefficient between particles and walls
CHINA PARTICUOLOGY Vol. 3, No. 6, 2005
362 μp
sliding friction coefficient between particles
μ ij
sliding friction coefficient between particles i and j
μr,ij
rolling friction coefficient, m
μ
rotational stiffness, m·N·s·rad
' r,ij
-1 -1
ωi
angular velocity of particle i, rad·s
ωnij
component of the vector of the relative angular velocity -1 of particles i and j in their contact plane, rad·s
Ω
given by Ω = ( r , t ) r ∈ Ω p , t ∈ T
Ωp
domain in which quantities are averaged, given by
{
{
Ω p = r r ≤ Lp
}
}
∂Ωp
boundary of domain Ω p
νi
Poisson’s ratio of particle i
σ1
principal stress, Pa
σn τ Π
compressive normal stress, Pa shear stress, Pa
∂Π
domain in R 2 boundary of Π
∂Π a
boundary with rigid constraints
∂Π t
boundary with prescribed traction
∂Π f
boundary with frictional force
φ
angle of material internal friction Dirac delta function
δ
References Alonso-Marroquin, F. & Herrmann, H. J. (2005). Investigation of the incremental response of soils using a discrete element model. J. Engng. Math., 52, 11-34. Babic, M. (1997). Average balance equations for granular materials. Int. J. Eng. Sci., 35, 523-548. Blight, G. E. (1986). Pressures exerted by materials stored in silos: Part I coarse materials; Part II fine powders. Geotechnique, 36, 33-56. Campbell, C. S. (1990). Rapid granular flows. Annu. Rev. Fluid Mech., 22, 57-92. Collins, J. (1990). Plane strain characteristics theory for soils and granular materials with density dependent yield criteria. J. Mech. Phys. Solids, 38, 1-25. Cox, A. D., Eason, G. & Hopkins, H. G. (1961). Axially symmetric plastic deformation in soils. Philos. Trans. R. Soc. Lond., A, 254, 1-45. Cundall, P. A. & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29, 47-65. de Gennes, P. G. (1999). Granular matter: a tentative view. Rev. Mod. Phys., 71, S374-382. Drescher, A. & de Josselin de Jong, G. (1972). Photoelastic verification of a mechanical model for the flow of a granular material. J. Mech. Phys. Solids, 20, 337-351. Farley, R. & Valentin, F. H. H. (1968). Effect of particle size upon the strength of powders. Powder Technol., 1, 344-354. Geng, J. F., Longhi, E., Behringer, R. P. & Howell, D. W. (2001). Memory in two-dimensional heap experiments. Phys. Rev. E, 64, 060301. Goldshtein, A. & Shapiro, M. (1995). Mechanics of collisional motion of granular materials, Part 1, general hydrodynamic
equations. J. Fluid Mech., 282, 75-114. Herrmann, H. J. & Luding, S. (1998). Modeling granular media on the computer. Continuum Mech. Therm., 10, 189-231. Hill, J. M. & Wu, Y. H. (1991). Kinematically determined axially symmetric plastic flows of metals and granular materials. Q. J. Mech. Appl. Math., 44, 451-469. Hill, J. M. & Wu, Y. H. (1992). Some axially symmetric flows of Mohr-Coulomb compressible granular materials. Proc. R. Soc. Lond., A, 438, 67-93. Hill, J. M. & Wu, Y. H. (1993a). Plastic flows of granular materials of shear index n: Part I yield functions. J. Mech. Phys. Solids, 41, 77-94. Hill, J. M. & Wu, Y. H. (1993b). Plastic flows of granular materials of shear index n: Part II plane and axially symmetric problems for n=2. J. Mech. Phys. Solids, 41, 95-115. Hill, J. M. & Wu, Y. H. (1996). The punch problem for shear index granular materials. Q. J. Mech. Appl. Math., 49, 81-105. Hill, J. M., Shi, J., Tordesillas, A. & Wu, Y. H. (1999). The velocity field for the punch problem for dilatant granular materials. Q. J. Mech. Appl. Math., 52, 99-110. Kishino, Y. (Ed.) (2001). Powders and Grains (a conference proceeding containing numerous DEM studies on the physics of granular matter). Lisse (The Netherland): A. A. Balkema Publishers. Lade, P. V. (1977). Elasto-plastic stress strain theory for cohesionless soil with curved yield surface. Int. J. Solids Struct., 13, 1019-1035. Langston, P. A., Tüzün, U. & Heyes, D. M. (1995a). Discrete element simulation of internal stress and flow fields in funnel 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: dependence of discharge rate and wall stress on particle interactions. Chem. Eng. Sci., 50, 967-987. Løvoll, G., Måløy, K. J. & Flekkøy, E. G. (1999). Force measurements on static granular materials. Phys. Rev. E, 60, 5872-5878. Lun, C. K. K., Savage, S. B., Jeffrey, D. J. & Chepurniy, N. (1984). Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J. Fluid Mech., 140, 223-256. Matsuoka, H. (1984). Deformation and strength of granular materials based on the theory of compounded mobilized planes and spatial mobilized plane. In Shahinpoor, M. (Ed.), Advances in the Mechanics and the Flow of Granular Materials (2, 814-836). Houston (USA): Bulf Publishing Company. Nedderman, R. M. (1992). Statics and Kinematics of Granular Materials. New York: Cambridge University Press. Nedderman, R. M. (1995). The use of the kinematic model to predict the development of the stagnant zone boundary in the batch discharge of a bunker. Chem. Eng. Sci., 50, 959-965. Oda, M. & Iwashita, K. (1999). Mechanics of Granular Materials. Rotterdam (The Netherland): A. A. Balkema Publishers. Oda, M. & Iwashita, K. (2000). Study on couple stress and shear band development in granular media based on numerical simulation analyses. Int. J. Engng. Sci., 38, 1713-1740. Polderman, H. G., Boom, J., de Hilster, E. & Scott, A. M. (1987). Solids flow velocity profiles in mass flow hoppers. Chem. Eng. Sci., 42, 737-744. Potapov, A. V. & Campbell, C. S. (1996). Computer simulation of hopper flow. Phys. Fluids, 8, 2884-2894. Ristow, G. H. (1992). Simulating granular flow with molecular dynamics. J. Phys. I, 2, 649-662.
Zhu, Wu & Yu: Discrete and Continuum Modelling of Granular Flow Rombach, G. & Eibl, J. (1995). Granular flow of materials in silos. Int. J. Bulk Solid Handling, 15, 65-70. Roscoe, K. H., Schofield, A. N. & Wroth, C. P. (1958). On the yielding of soils, Geotechnique, 9, 71-83. Rothenburg, L. & Selvadurai, A. P. S. (1981). A micro-mechanical definition of the Cauchy stress tensor for particulate media. In Selvadurai, A. P. S. (Ed.), Proceedings of the International Symposium on the Mechanical Behavior of Structured Media (pp.469-486). Ottawa, Canada. Schmidt, L. C. & Wu, Y. H. (1989). Prediction of dynamic wall pressures on silos. Int. J. Bulk Solid Handling, 9, 333-338. Seville, J. P. K., Tüzün, U. & Clift, R. (1997). Processing of Particulate Solids. Landon: Blackie Academic & Professional. Shen, H. & Hopkins, M. (1988). Stress in a rapid, simple shear flow of granular materials with multiple grain sizes. Part. Sci. Technol., 6, 1-15. Spencer, A. J. M. (1964). A theory of the kinematics of ideal soils under plane strain conditions. J. Mech. Phys. Solids, 12, 337-351. Spencer, A. J. M. (1982). Deformation of ideal granular materials. In Hopkins, H. G. & Sewell, M. J. (Eds.), Mechanics of Solids (pp. 607-652). New York: Pergamon Press. Steingart, D. A. & Evans, J. W. (2005). Measurements of granular flows in two-dimensional hoppers by particle image velocimetry. Part 1: experimental method and results. Chem. Eng. Sci., 60, 1043-1051. Stewart, R. L., Bridgwater, J., Zhou, Y. C. & Yu, A. B. (2001). Simulated and measured flow of granules in a bladed mixer ⎯ a detailed comparison. Chem. Eng. Sci., 56, 5457-5471. Tanaka, K., Nishida, M., Kunimochi, T. & Takagi, T. (2002). Discrete element simulation and experiment for dynamic response of two-dismensional granular matter to the impact of a spherical projectile. Powder Technol., 124, 160-173. Walton, O. R. & Braun, R. L. (1986). Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheol., 30, 949-980. Wightman, C., Moakher, M., Muzzio, F. J. & Walton, O. R. (1998).
363
Simulation of flow and mixing of particles in a rotating and rocking cylinder. AIChE J., 44, 1266-1276. Wu, Y. H. (1990). Static and Dynamic Analyses of the Flow of Bulk Materials Through Silos. PhD thesis, The University of Wollongong, Australia. Wu, Y. H. & Schmidt, L. C. (1992). A boundary element method for prediction of silo pressure. Int. J. Comput. Struct., 45, 315-323. Xie, H. Y. & Shinohara, K. (1999). Measurement of solids velocity in a conical hopper by mass tracer particles. Chem. Eng. Sci., 54, 455-459. Yamane, K., Nakagawa, M., Altobelli, S. A., Tanaka, T. & Tsuji, Y. (1998). Steady particulate flows in a horizontal rotating cylinder. Phys. Fluids, 10, 1419-1427. Yang, R. Y., Zou, R. P. & Yu, A. B. (2003). Microdynamic analysis of particle flow in a horizontal rotating drum. Powder Technol., 130, 138-146. Zhang, Y. & Campbell, C. S. (1992). The interface between fluid-like and solid-like behaviour in two-dimensional granular flows. J. Fluid Mech., 237, 541-568. Zhu, H. P. & Yu, A. B. (2001). A weighting function in the averaging method of granular materials and its application to hopper flow. Bulk Solids Handling, 21, 53-57. Zhu, H. P. & Yu, A. B. (2002). Averaging method of granular materials. Phys. Rev. E, 66, 021302. Zhu, H. P. & Yu, A. B. (2003). The effects of wall and rolling resistance on the couple stress of granular materials in vertical flow. Physica A, 325, 347-360. Zhu, H. P. & Yu, A. B. (2004). Steady-state granular flow in a 3D cylindrical hopper with flat bottom using DEM simulation: microscopic analysis. J. Phys. D: Appl. Phys., 37, 1497-1508. Zhu, H. P. & Yu, A. B. (2005a) Micromechanic analysis of hopper flow of granular materials. J. Engng. Math., 52, 307-320. Zhu, H. P. & Yu, A. B. (2005b). Steady-state granular flow in a 3D cylindrical hopper with flat bottom using DEM simulation: macroscopic analysis. Granular Matter, 7, 97-107. Manuscript received November 2, 2005 and accepted November 10, 2005.