Efficient computer simulation of moving granular particles

Powder

Technology,

78 (1994)

Efficient computer Caroline Depariment

(Received

Hogue

51

51-66

simulation of moving granular

and David

of Engineering,

Universiv

particles

Newland of Cambridge,

Trumpington

Street,

Cambridge

CB2 IPZ

(UK)

March 1, 1993; in revised form June 25, 1993)

Abstract A two-dimensional model for computing contacts and motions of granular particles of different shapes, sizes and material properties is presented. The primary aim of this model is to achieve a high degree of computational efficiency, to allow simulations to be performed very rapidly on a modest sequential machine. The important features of the model are (i) a polar representation of particle shape, (ii) a simple algorithm for contact detection, (iii) a binary collision model using momentum exchange as the describing feature, (iv) energetic coefficients of restitution selected pseudo-randomly, (v) a Coulomb friction approximation to describe slip at impact, and (vi) an optimized move-update algorithm. Results in satisfactory agreement with analytical solutions and existing, less efficient computer simulations are presented.

1. Introduction

Although several models describing the motion of particles in two dimensions have already been formulated, none of them provides sufficient computational efficiency to run large simulations in reasonable time scales. We define reasonableness according to a time-benefit trade-off. We feel this excludes all approaches which require more than about 2 h of computing time or specialized hardware. This contrasts with the time scales of 8 h on a Cray supercomputer [l], or even days of computation with previous models [2], most of which (including the Monte Carlo method [3]) deal only with spherical particles. This article discusses a new approach for handling particles of arbitrary shape. The proposed momentumexchange model, which we call the contucf model, is described in detail in Section 3. The main areas of application of the contact model include vibrated screening and sieving, and hopper and chute flow. This model represents the flow at the microscopic (or local) level. Although many large-scale applications have been simulated using continuum models (fluid-like behaviour), these models use equivalent fluid properties derived experimentally and they cannot represent local effects, for example bridging in hopper flow. The contact model (small-scale) together with the partition model (global, large-scale) allows large-scale simulations to be performed on reasonable time scales while still showing localised behaviour patterns. The partition model is the subject of an earlier publication [4] and will not be discussed again here.

0032-5910/94/$07.00 0 1994 Elsevier 0032-5910(93)02748-Y

SSDI

Sequoia.

All rights reserved

In this paper, existing models for simulating the motion of granular media are reviewed briefly in Section 2. Our contact model and its limitations are then described in detail in Section 3. Typical results, which include rotating blocks and packing simulations, are given and discussed in Section 4. The paper concludes with a discussion of further ideas for research in the computer modelling of granular motion.

2. Existing models The existing approaches to modelling granular materials can be divided in two categories: continuum and discrete models. Continuum models assume the granular media to behave as a fluid. There is extensive literature on continuum approaches (e.g., refs. 5-7). These approaches are mathematically complex, and yield results that sometimes differ by up to an order of magnitude from experimental values [8]. Currently continuum approaches are limited to the modelling of spherical granular particles. Discrete models, on the other hand, are conceptually much simpler but rely heavily on computing power when simulating medium- to large-scale systems. The discrete approaches can be divided into three main classes of model: (i) classical Newtonian dynamics models; (ii) statistical mechanics models; and (iii) hybrid models. These categories are discussed briefly in the following sections.

52

2.1. Classical Newtonian dynamics approach

This approach involves the application of Newtonian dynamical equations to a system of impacting particles. This includes keeping track of all the forces and moments acting on each particle at every time step, and integrating the equations of motion to obtain the new state of the system at the end of each step. We have tested a system of one hundred spherical particles [2], where the contact forces were computed assuming a spring and dashpot to model the interaction of individual particles (in both normal and tangential directions). The results obtained agreed satisfactorily with experiments [2]. However, these simulations were achieved only at great computational cost. The time scales required were well outside the desired ranges; to produce the simulations shown in ref. 2, a few weeks of computation were necessary on an HP-800 workstation. Cundall and Strack’s discrete element model (DEM) [9] is probably the first and best-known formalized model based on the computation of forces and the solution of Newtonian dynamical equations for nonspherical particles. Successful results were published for soil and rock mechanics applications using polygonal particles (or polyhedral blocks in three dimensions) [lo, 111. Other examples of applications based on the DEM include the study of granular flow by Walton [12], vibrationally induced segregation by Haff and Werner [13] and plug flow by Tsuji et al. [14]. As a research tool, where time constraints are not crucial, DEM approaches have proved valuable. Again, however, the large computational resources required appear to prevent the method’s widespread use as a design tool or in industrial applications. Another important limitation of DEM approaches is the difficulty of handling both spherical and polygonal particles in the same simulation. This difficulty is due to the different particle representations and contact detection algorithms used for spherical and polygonal particles.

2.2. Statistical mechanics approach An alternative to computing all the forces and integrating the equations of motion for an entire system is to use a statistical approach. This eliminates the need for force calculations and integration routines, as the state of the particles is determined according to an assumed statistical distribution. Two such statistical models are based on the Monte Carlo ideas of statistical mechanics [3,15]. These models are simple and efficient. Unfortunately, they are also inflexible: there is little control over most parameters, e.g., the particles can be modelled only as hard spheres, and the vibrations - if present - must be harmonic. Also, these models tend to apply only to a specific application.

2.3. Hybrid approach Some models adopt a hybrid approach, combining ideas from both deterministic and statistical approaches. These include Hawkins’s approach [ 161,which simulated the granular flow of 2-D rigid disks by using stochastic boundaries and binary collisions. The simulations were relatively efficient, and they were performed on an advanced personal computer. Another example is the Goldreich-Tremaine model [17] for the study of the motion of particles in planetary rings. The occurrence of a contact is determined statistically according to a Boltzmann distribution, but the outcome of a contact is computed using energy and momentum considerations. A more general hybrid model, applicable to particles of arbitrary shape, is the contact model described in the next section.

3. A new model: the contact model A simplified algorithm of the model is depicted by the flow chart in Fig. 1. The two main stages of computation are (i) the contact processing and (ii) the motion update. In the contact processing stage, the state of the particles after impact is computed using only energy considerations and momentum conservation laws. In order to do this efficiently, all collisions are assumed to be binary. That is, particles do not impact with more than one neighbouring particle simultaneously. Thus, if the contact detection routine detects a multiple contact, it is not interpreted as such, but as a series of successive impacts all occurring rapidly within the given time step. If this assumption of binary collisions

Process

contacts

Check trial configuration for overlaps and reorder particles

Fig. 1. Simplified

flow chart for the contact

model.

53

is not enforced, the equations of conservation of momentum become much more complex, decreasing substantially the efficiency of the computation. Using the assumption of practically instantaneous but sequential contacts, particles can be considered to be in free flight between impacts. This is the main assumption of the motion update stage, and avoids the need for integration. However, free-flight motion is not possible in all cases. For example, free flight cannot occur when a particle is surrounded by others in such a way that a displacement in the direction of its resulting velocity generates an excessive overlap with a neighbouring particle. In such a case, the particle may not be moved. To overcome this problem, the free-flight motions of the particles are considered only as trial configurations (similar to the trial configurations of the Monte Carlo method [3]). If the trial state does not involve an excessive overlap, it is accepted; otherwise the old position is retained as the new position. The rate of acceptance of the trial configurations, using the method described in Section 3.7, is optimized to ensure a minimum number of rejected configurations per time step (see Section 3.7.2). 3.1. Assumptions As can be seen from the above brief description of the contact model, a number of simplifying assumptions have been made. This is because of the large number of factors that influence particulate behaviour in systems of non-ideal particles. Simplifying assumptions have been found to improve greatly the computational efficiency without significantly affecting the overall particle behaviour; hence the utility of the model. In detail, the assumptions are as follows: (i) any particle shape is represented by a series of vertices joined by straight line segments; the restrictions of the possible shapes are explained in Section 3.3; (ii) each granular particle retains its shape after impact (quasi-rigid particles); (iii) motion is two-dimensional with the centre of mass of all particles moving in one plane; (iv) each impact is instantaneous and binary; if more than one contact is detected involving the same particle, the contacts are assumed to happen sequentially, with the contact having the highest relative velocity taken as the first contact; (v) contacts occur at a point; this allows the angular momentum of both particles individually to be conserved about the contact point; (vi) pseudo-random energetic coefficients of restitution (which vary at each impact according to a normalized random distribution centred at the average value) are used in the energy dissipation equations in the normal direction of impact;

(vii) Coulomb’s law of friction determines the energy dissipation in the tangential direction of impact; and (viii) Stronge’s energy dissipation hypothesis is assumed as the basis for determining energy dissipation in the normal direction; this is described in detail in Section 3.6.1. 3.2. Critical time step For models with elastic impacts, where particles have stiffness and damping, the minimal time step required to ensure numerical stability is related to the minimum contact duration. That is, several integration steps must be performed for each contact. In the literature, the minimum contact duration varies from 3 time steps [18] to 10 time steps [12]. Most discrete element applications use an approach similar to Cundall’s, where time step of the critical the determination (tcrit =2&Z) is based on a single-degree-of-freedom system of mass m connected to the ground by a spring of stiffness k. A factor of safety is normally applied to this value, so that the actual time step is always less than the critical value. For models with instantaneous impacts, such as the contact model, the minimal time step is usually based on the minimum time between inter-particle collisions. This minimum time is thus related to the maximum relative velocities between particles during a simulation. A suitable approximation for the maximum relative inter-particle velocity, and thus a limiting time step, can be obtained for a given simulation as follows. Assume that the particles are initially at rest. If there are no external forces applied to the system other than gravity, the maximum velocity of a particle is proportional to the maximum drop in height H V,,, = L&z

(1)

so that the maximum relative inter-particle velocity (Av,,,~) can be approximated by twice this value: Av,,, = 2m

(2)

The maximum displacement of a particle having the above velocity during a time step can be written as A~max= -

IgAt2+2mAt n

L

The limiting value that we impose on this maximum displacement is that it must be less than the minimum radius (in the sense of our polar representation) to ensure that the maximum overlap detected is never more than half the smallest diameter of a granular particle. Neglecting higher orders in At, we thus get A~max
or, for At -

(4) 0,

54

2cH

At < R,i,

(9

which yields a critical value for the time step of Atcrit= -

Rmin

(6)

2433

As an example, if the maximum vertical drop is 1 m, and the minimum radius is 10 mm, the critical time step is then Att,,i,=O.OO1s. In the case where there is an external force applied to the system of particles, a vibrating floor for example, the maximum energy contribution of the external force is included in the calculation of the maximum interparticle velocity, and the derivation of the minimum time step using this new velocity is performed as above. 3.3. Representing arbitrary particles One of the advantages of the contact model over many existing models is that it can handle particles of different shapes. Particle shape is represented by vertices joined by small line segments, as shown in Fig. 2. This information is stored in a 2-D array in polar form. That is, the ‘radius’ from the centre of mass of each vertex is stored along with its angle computed from a reference point at 0 radians. The number of vertices used varies for each type of particle; an average of about 24 vertices is used for the particles in the simulations presented in this paper. A polar representation helps in the contact detection process, as will be discussed later on. Note that this representation also has the advantage of not restricting the shapes to be convex, as is the case with some particle representations. Although some shapes cannot be represented by the above scheme, for example a narrow particle in the shape of the letter ‘c’, the representation of all the common granular material shapes is straightforward. In addition to shape information, the data structure of each particle contains information on density, moment of inertia about its centre of mass, averaged coefficient of restitution, and friction coefficient (between particles and walls), as well as motion information, i.e., positions and velocities.

__--_-

Fig. 2. A particle

3.4.1. Inter-particle con tact detection The detection of inter-particle contacts is performed according to a simple geometric algorithm. At present, all particles are checked for possible contacts. A technique assessing only neighbouring particles for possible contacts is described in ref. 19 and will be implemented in future versions of the model. This technique produces a list of neighbours for each particle, and looks only for contacts between these neighbours. The list of neighbours is not updated at each time step; it is updated only when a particle’s position has changed by more than a given amount. Thus, for high-energy systems, where particles are moving at high velocities, the technique in ref. [19] may not provide much increase in efficiency compared with checking all particles for contacts. The efficiency increase is more significant in dense, slow-moving systems. The algorithm for contact detection, as implemented, is given in Fig. 3 and works as follows. The first stage consists of identifying possible contact pairs. Since the maximum distance (R,,,) of each particle to its centre of mass is stored in its data structure, the simple check (R,,,i + R,,J

> cm-distance

(7)

where cm-distance is the distance between the centres of mass of particle i and j, identifies a possible contact. For systems consisting only of spherical particles, this condition is sufficient to identify a contact. Notice that this check quickly rejects many potential contact possibilities. The next step is a more detailed analysis to determine the contacts precisely. This relatively costly analysis is performed only for particles for which condition (7) is satisfied. The search is done as follows. Each vertex of particle j, which belongs to the set defined by the shaded region shown in Fig. 4, is checked to see if it lies inside particle i, until an overlap is detected. This for all pairs of neighbouring particles i and j Identify possible collisions by comparing relative positions of centres of mass (cm-distance) with the sum of maximum radii. if collision possible then Perform detailed contact check: each vertex VEI of particle i is checked to see if it lies inside particle j and vice versa, until greatest overlap is found. end for

Oradian

represented

3.4. Contact detection The detection of contacts between a particle and a wall is done in a straightforward fashion, where the vertices of each particle in the vicinity of a wall are checked for overlap with a wall. The detection of interparticle contacts is discussed below.

in

polar form

using 16 vertices.

Fig. 3. Algorithm

for interparticle

contact

detection.

55

points in this region need to be checked

Only peripheral

Fig. 4. I is the set of all vertices

lying in the shaded

region.

i

1

(4 Fig. 5. (a) Typical contact between of particle i containing point p.

(b) particles i and j; (b) triangle

process is repeated for each vertex of particle i for possible overlap with particle j. This second check is done to identify the deepest (largest) contact. Figure 5(a) illustrates a typical contact. Let d be the distance between vertex p of particle j and the centre of mass of particle i. From the angle that the line d makes with the horizontal, the triangle of particle i in which p may lie can easily be found. Thus the values of R, and R,+l as shown in Fig. 5(b) can be determined from the data structure. A simple check to determine whether p lies inside the triangle is to compute areas. If the area of triangle iAB is greater than the sum of the areas of the two triangles iAp and ipB, then p lies inside particle i. In the event that more than one vertex of particle j lies inside particle i, or vice versa, the one producing the deepest contact that is moving towards the other particle is retained. In the event that the velocity at the points of contact is not in the direction of the other particle for any vertex, then the deepest contact is retained. The depth of a contact is given by the distance between point p and the edge particle i in the direction normal to the impact. The determination of this normal is described in Section 3.5. A discussion on what is meant by particles moving towards each other (or away from each other) at the contact point is given in Section 3.6.

3.5. Normal direction of impact Finding accurately the tangential and normal directions of impact, represented by the X-Y coordinate system in Fig. 6, is crucial to the contact processing stage, as it was found by experimentation that errors of more than 10 degrees may lead to non-realistic particle motion. In our particle representation, the determination of the normal is a geometric problem. When a vertex penetrates a flat surface, the normal direction of impact is always assumed to be perpendicular to the flat surface (Fig. 7(a)). When, on the other hand, a vertex penetrates a comer, as shown in Fig. 7(b), the normal is chosen perpendicular to the line joining the vertices on either side of the corner (segment m-n in Fig. 7(b)). In order to differentiate accurately between the above two types of contacts, the following algorithm is used. Let p be the vertex chosen as the contact point on particle i. If p is the only vertex lying inside particle j, the situations are comparable to those shown in Figs. 7(a) and (b). Then, by checking which segments (i.e., which of wrn, m-o, o-n or n-w) are crossed by the lines 4-p and p-r, it is possible to identify the type of contact correctly. In the case where more than one vertexp of particle i lies inside particle j (see Fig. 7(c)), the points 4 and r are moved to the first two vertices that lie outside particle j, on either side of point p. The normal is again taken to be perpendicular to the line joining m and n. As for the case where particle i overlaps with more than one point of particle j, the two points farthest apart of particle j that lie inside particle i are taken to be m and n, and the normal is taken to be perpendicular to the line joining them, as illustrated in Fig. 7(d). Finding the correct type of contact is rapid and requires little computational time, but is important for the contact processing stage.

ty

isi? I

j (XilYi) I@

_--_- px

_____

PY

~wz

-+

X

I

Fig. 6. The impulse components acting on body i in the normal and tangential directions of impact.

56

\

i

i

i

’

n9 1’0

rn\,

V

r

\

Fig. 7. Normal direction of impact: (a) for one-point contacts onto flat surfaces; multiple-point contact of i inside j; (d) for multiple-point contact of j inside i.

3.6. Contact ordering and update

When all the contacts have been identified, an update of the linear and angular velocities of each particle involved in a contact is performed according to a specific ordering. This ordering applies to all contacts. The contacts are ordered according to the closing velocity at the contact point, highest velocity contacts first. The ‘velocity of closing’ of a contact is given by the difference in velocity in the normal direction, which is oriented such that a positivevelocity difference implies that the two particles are moving towards each other. If, for example, particle i is involved in multiple contacts, the closing velocity at each contact point is computed. The fastest-moving contact is then processed. This processing yields a new velocity for particle i, and so a new check of the closing velocities is performed. Again, the fastest moving contact is updated, and so on. Since this is done sequentially, not all contacts may be seen as moving towards each other by the time they are processed. Only those moving towards each other require processing. 3.6.1. Interparticle contact update The contact model ensures that linear and angular momentum are conserved during each impact. Energy is dissipated in the normal direction of impact through a pseudo-random coefficient of restitution and in the tangential direction through friction. Using the X-Y coordinate system shown in Fig. 6, the equations of linear momentum can be written as mi( Ye- nix>= Px

(8)

Illi( l& - Q) = P,

(9)

rnj( l& - Q) = - Px

(10)

mj(z$-v,~)=

(11)

-P,

where Px and PY are the X and Y components of the impulse at impact, v, and vjy are the components of the velocity of the centre of mass of body i before impact, and v& and viy are the corresponding velocities after impact.

r

P n\. %

PW (b) for one-point

W

i

%

(c)

n

m

\ VJ

!n

I I I

I

\ V

4

(4 contacts

P onto corners;

(c) for

Since contacts are assumed to occur at a point, conservation of angular momentum about the point of impact gives (ri--r,,)mi~~+ziwi=(ri-r,,)m,v:+Iiw:

(12)

(~j-r,,)mi~+ljwj=(rj--r,,)mjV~+ljo~

(13)

where ri and r,, are the position vectors of the centre of mass of body i and the point of contact, respectively, mi and Ii are the mass and moment of inertia of body i and w, and 6~: are its angular velocities before and after impact. These angular momentum equations can also be written in terms of the impulse components as Ii (wi - CL+) = Pxyi - Pyq

(14)

h(&.$- Wj)= - Pxyi + PYXj

(15)

where xi and vi are the distances from the centre of mass of body i to the contact point along the tangential and normal directions of impact. Several approaches were considered in order to model the energy dissipated in the normal and tangential directions of impact. These are described in detail in ref. 20. The approach adopted is given below.

Inter-particle friction and restitution. The modelling of the energy dissipated through restitution and friction at impact is done using Stronge’s energy dissipation hypothesis 1211to govern the energy losses in the normal direction, whilst Coulomb’s friction law models the energy losses in the tangential direction. Stronge’s energy dissipation hypothesis explicitly relates the coefficient of restitution (called the energetic coefficient e*) to the energy dissipation by an irreversible deformation process. This hypothesis results in non-frictional dissipation being proportional to (1 - e’, ). Thus, the total impulse is separated into parts for each phase of unidirectional slip, and the response is computed using only two independent material constants, e, and p. This hypothesis has been adopted as it avoids the spurious energy gains during impact

which may occur when using Newton’s law of restitution or Poisson’s hypothesis* (see ref. 21). The energetic coefficient is defined as follows. The square of the energetic coefficient is the ratio of elastic strain energy released at the contact point during restitution to the energy absorbed by internal deformation during compression: 2_ e*--

I&- IV, _ WY

(16)

where I?‘,, is the work done by the normal component of impulse during compression, and I$JYis the total work done by the normal impulse. This definition is equivalent to both Newton’s kinematic coefficient of restitution and Poisson’s kinetic coefficient of restitution unless the bodies are rough, the contact is eccentric and the initial sliding velocity halts or reverses during the impact. In order to use this energetic coefficient, analysis of the slip modes at each impact is required. Analysing the phases of slip. A useful method for analyzing the different phases of slip is due to Routh (see ref. 22). In Routh’s method, the changes in the velocities during a collision are given as a function of the normal and tangential components of the impulse, P, and Pz The collision process is divided into a period of compression (bodies approaching each other) and a period of restitution (bodies moving apart). For convenience, the notation used here follows that of Wang and Mason [22]. For a contact between bodies i and j, the components of the impulse Px and PY, shown in Fig. 6, are given by eqns. (8) to (11). The tangential and normal components of the relative velocity at the point of contact are denoted by S (for sliding) and C (for compression) and are given by

dC =sjd3+B2 dp,

where s =sgn(S,) = ]&]/So, and B1, B, and B, depend on the inertia of the system and are given by B1=l,L+$+f mi

(21) q

,

J

B,‘+I+yf+y: 2

mi

mj

Ii

(22)

4

(23) From eqns. (19) and (20), we can obtain the range of parameter values for a closing contact (where the relative normal velocity is such that the two bodies are moving towards each other): -&
(24) Integrating eqns. (19) and (20) for the initial period of sliding gives the following conditions for So # 0: (i) for sliding to continue throughout the impact: ,. (B2+sd3,WO

>

(s&

+

&I(

-

CO)

2

(25)

(ii) for sliding to halt or reverse duiing compression: (& +s/&)s& < (s& + 4.4 )( - G) (26) and (iii) for sliding to halt or reverse during restitution: (& + @I)( - G) < (S +sILB&& ,. < (s&+&)(-G)

2 Y

(27)

where py is the terminal normal impulse at separation and P, is the normal impulse at the transition from compression to restitution. That is, St Py=

s

F(t) dt

(28)

0

where v,, and v,,, are the components of the velocity of body i at the contact point. Note that S,, and C,, are the initial values of the sliding and compression velocities. The changes in the normal and tangential components of the impulse during sliding are related by Coulomb’s law of friction, dPx= -ru. dP, sgn(S). Using Coulomb’s law, the equations of relative motion of the bodies can be written as *In Poisson’s hypothesis, the coefficient of restitution is defined as the ratio of the normal impulse during restitution to the normal impulse during compression.

with py and py corresponding to the two different values of the upper limit 6, at the end of the impact and at the end of compression. Also, if sliding halts, reversed sliding will occur (both in compression and restitution) provided that /_L
For a direct impact (S,=O), only two modes of slip can occur: continuous sliding and sticking during compression [22]. These are determined according to the value of the coefficient of friction as follows: p
(30)

58

I_L > II, for sticking

(31)

where pL, is as defined above. Analytical solutions of imp&e.

-co

I’,= B2

The equations

for the impulses using Stronge’s hypothesis are obtained using the same method as used by Wang and Mason [22] (see the Appendix for a detailed derivation of these equations). The collective equations for the terminal impulses for the different slip modes using the energetic coefficient approach, and where s=sgn(S,), are (i) sliding: &=

where

-s&,

+ wB3

c0 &+V&

B3fs1.4

p =p _ Y Y

-s@3)

[ py-P-y,-(

1 &+vJ% 1

py$;

(37)

[

(B1_s$$iftsti,j2

1

C + &P&C0

E+

1 l/2

(38)

where B,--s@,

[

’

B3+s/-4

1

-C”

I’,=

l/2

B2

(35)

2

(44)

where

(36)

Py=Py+e*

(d-F;)2-~*P$~)1~] 3

(34)

‘“,+

-

py-Py)2

2s

2&

rj,-

py=py+e*

p,=

(B2+sd3) 032

-

-

”

(43) )

(33)

[

”

+v&

2so

Ij,= --s/J P,-

(ii) sticking in compression: 2

(42) B3

(32)

-(l+e*)

&=

&

P,=

(v) reversal in restitution:

x

F,=

and

(39)

(iv) sticking in restitution: (40)

and

+ ~43

Py=

so B3+wh

(45)

Contact update. The velocities of each body at the end of the impact are obtained by substituting the values of the impulse components in eqns. (8) through (11) and in the equations for the conservation of angular momentum at the point of contact, eqns. (14) and (15). Thus there is no need to solve a system of equations. It should be noted that the value of the energetic coefficient of restitution e, for each body is chosen pseudo-randomly according to a Gaussian distribution centered at an average experimental value C?* , with an experimental standard deviation u. This randomness is included to account for the fact that e, is not constant but is highly dependent on geometry and impact velocity as well as on materials. In addition, a composite coefficient of restitution [23] is computed for contacts between bodies i and j having respective coefficients of restitution e *i and e *i, and is given by

e’, =

mje2,

mi+mj

i

+

mi&.

j

(46)

m,+T?lj

3.6.2. Floor and wall contact update

-

(41)

For the contacts between a particle and a wall or floor, the same equations apply. The mass of a wall or floor is considered to be infinite; hence its velocity is unaffected by impacts. Once the slip mode has been identified using inequalities (25)-(27), the impulse com-

59

ponents at impact are computed using the appropriate equations given in the previous section. Equations (8), (9) and (14) are then used to obtain the final velocities directly. 3.7. Particle motion At the end of each cycle of computation (time step), a trial configuration of new particle positions is generated. It is based on the assumption that the particles are in free flight (with gravity as the only external force) between each impact. 3.7.1. Trial co&uration A trial configuration is obtained as follows. Since the state of each particle after its last impact is known, and since the free-flight assumption implies that gravity is the only force acting on the particles, explicit equations can be derived for the displacement and velocities of each individual particle. So for each particle i, the displacement equations can be written as xi = g At2 sin p +ii At +xj 2

(47)

y:= - ; At2 cos p +ji At +yi

(48)

O:=a At+@

(49)

and the velocity equations

as

$ =g At sin p +&

(50)

r’:= -g At cos p+$

(51)

8; = fji

(52)

where xi, yi and Oiare the X-Y and angular components of the displacement of the centre of mass of body i; &, ji and ii are the X-Y and angular components of its velocity; and p is the angle between the X axis and the horizontal. The prime indicates the trial (updated) values, and At is the size of the time step. These explicit solutions provide a quick way of obtaining a trial configuration for each particle of the system. In order to capitalize on this efficiency, the evaluation of this trial configuration must also be efficient. The method used for the evaluation and the determination of the new state is described below. 3.7.2 Motion update This assumption of free flight generates only trial configurations, since, for example, a particle may be in a position where any movement in the direction of the resulting velocity, or any amount of rotation, would create an overlap with one or more surrounding particles, or an excessive overlap with the floor or walls. In such cases, the old positions are retained as the

new ones, but it should be noted that the velocities are always updated at the contact processing stage if the particles are involved in a contact. Although the (small) movement that would otherwise occur is neglected, the change in velocity (and therefore momentum) is represented correctly. This ensures the best accuracy within the limitations of the model. The process of determining which trial configuration must be rejected involves a partial ordering scheme that tries to optimize the acceptance rate of the trial configuration without compromising the overall efficiency. An exhaustive check of all permutations of particles’ current and trial positions (a total of (2”)! permutations) would be intractable. Instead, a heuristic algorithm is used to try to maximize the number of particles moved. The algorithm used is summarized in Fig. 8. The first step consists of finding all the contacts that occur if all particles are moved to their trial configuration. For every such contact found, at least one particle position must later be rejected if this overlap is to be avoided. The second step detects contacts between each new trial configuration and the current position of all other particles. This information allows us to fill in a table of partial orderings, keeping track of which particle must be moved before another to avoid an overlap. Figure 9 illustrates this procedure with a simple example. In Figs. (9a)-(c), the solid lines represent the current positions A, B, C and D (before update), and the broken lines the trial positions A’, B’, C’ and D’. From Fig. 9(b), we can see that the trial configurations cannot all be accepted, since the contact B’-C’ exists. In step 4, to decide which particle(s) to reject so as to minimize the number of rejections, the contacts between each trial configuration and all other current positions are recorded, as shown in Fig. 9(c). In this example, the following information is found: B4A B
(53)

which can be written in matrix form as

1. Detect contacts for entire trial configuration 2. Detect contacts between each trial configuration and all other current configurations 3. Compute partial order < to prioritize particles 4. Apply Warshall’s algorithm to compute -c+, the transitive closure of + 5. For all contacts detected in (l), select a particle to reject (the one with least dependencies or one already rejected) 6. Reject all dependencies of particles rejected in (5) 7. Update contact list Fig. 8. Motion update

algorithm.

I

(cl

(4

Fig. 9. Example illustrating the steps of the motion update algorithm for particles A, B, C, and D: (a) current configurations (before update): (b) trial configurations; (c) each trial configuration vs. other current positions; (d) updated configurations.

(54) D 1 Thep,
A B C D

ABCD

the last column representing t:ae total number of dependent relations (or number of particles related to the one on that row) for each particle. Steps (5) and (6) consist of selecting the particles to be rejected. For each collision found in Fig. 9(b), at least one of the particles involved must be rejected. In order to reject a minimum number of particles, the particles with the lowest number of dependent relations are rejected. In our example, the only contact between trial configurations is B’-C’, and from the above table B has three dependent relations and C has none, so C is selected. If C had some dependent relations, these would also have to be rejected, to avoid any overlap. At this stage, a last check is performed to ensure that all rejected particles are recorded as being involved in a contact, so that their velocities are updated accordingly in the next contact update stage. All nonrejected particles keep their trial configuration as their new state, as shown in Fig. 9(d). The contact detection information gathered at this stage is used at the beginning of the next time step. Since, for interparticle contacts, no overlap is allowed when a contact is detected, all the required information about the contact is gathered as though the particles actually touch, even though their trial motion will be rejected. Moving froor overlaps. In simulations involving a moving floor, large overlaps may be created between particles and the moving floor. That is, although the force applied by the floor may be sufficient to move the particles upward, the nature of the model is such that the upward trajectory of the particles may lag slightly behind that of the floor, depending on possible inter-particle overlaps in the trial state. For example, if a particle already in contact with the floor has its trial state rejected, it will remain at its previous position even though the floor will continue to move. And if this movement is in the upward direction, it will create a theoretical overlap with the particle. In order to minimize such situations, a particle’s trial configuration is accepted even though an overlap may exist with the moving floor, as long as the particle is moving away from the floor. In this way, the motion-update strategy is applied in its simplest form.

4. Experimental

results using the contact model

no. dependents 1 3

(55)

In order to validate the contact model, analytical solutions of sample configurations and existing simulations for comparisons were sought. Applications ranging from falling dominoes to bin flow are discussed below.

61

4.1. Domino effect; steady propagation

In this set of simulations, the contact model is used to model the behaviour of the successive toppling of regularly spaced dominoes. This phenomenon, commonly called the domino effect, has been investigated by Stronge [26] and Stronge and Shu [27]. These two studies show that the wave of destabilizing collisions produced by the toppling of dominoes approaches asymptotically an intrinsic or natural speed, provided that the spacing between the dominoes is sufficiently large. Also, this intrinsic speed depends mostly on the spacing, the coefficient of restitution and the sliding friction during each collision. Experiments were reported and the results compared against both the single collision theory [26] and the cooperative group theory [27]. The cooperative group theory, which accounts for the effect of the toppled dominoes leaning against each other, proves more ‘accurate unless the dominoes are widely spaced, in which case the assumption of a single collision between neighbours is reasonable. Simulations using the contact model have been performed and compared against both these experimental and theoretical results. The parameters used in the simulations are the ones published by Stronge [26] and Stronge and Shu [27], and are given in Table 1. Figure 10 shows a schematic diagram of the dimensions and spacing of toppling dominoes. The simulations are performed using the three types of dominoes used by Stronge: plastic (OpaleneB) commercial dominoes, and Perspexa and TufnolQp blocks. In the simulations, the natural speeds were reached after the collision wavefront had passed through the first 4-15 blocks, depending on the spacing and the initial toppling speed. After this initial phase, the wavefront tended to settle down to a steady speed. A similar behaviour was observed by Stronge and Shu in their experiments. TABLE 1. Domino dimensions friction

and coefficients of restitution

Type of domino

L

h

(mm)

(mm)

Opalene Perspex Tufnol

41.8 80.0 80.0

7.5 9.9 9.6

Fig. 10. Toppling

Jt

1

2

3

4

5

.spdng, X/h

Fig. 11. Comparison of simulations with Stronge and Shu’s experiments and theory for natural speed of toppling dominoes. The types of domino used are as follows: Experiments: (0) Opalene dominoes; (0) Perspex blocks; (0) Tufnol blocks; Simulations: ( x ) Opalene dominoes; (+ ) Perspex blocks, ( *) Tufnol blocks. Theory: curve: h/L = 0.12, p = 0.2, e = 0.0.

The simulation results are compared with the predicted results of Stronge and Shu’s cooperative group theory and the experimental results shown in Fig. 11. The non-dimensional natural speed of the wavefront V*/@ for the dominoes in the simulations is in good agreement with the experimental values. Figure 12 shows snapshots of one of the simulations. High-speed photographs of a corresponding experiment reported by Stronge and Shu [27] are shown in Fig. 13. The contact model accurately shows that the interaction between every pair of neighbours is not limited to a single type; that is, some pairs of blocks at the wavefront have sliding contacts and others have single or multiple collisions. Agreement between simulations and experiments is best when N/I > 1.5. When the spacing is small (A/h< 1.5), the number of collisions between neighbouring blocks is large, so the chance of particle motions being rejected in the contact model is greater. Then the computed natural speed of the wavefront can be less than its experimental value. This can be observed in Fig. 11.

and

e

CL

0.85 0.55 0.65

0.18 0.25 0.15

of regularly spaced dominoes.

1.0

4.2. Pouring into a container Rosato et al. [3] have simulated the pouring of disklike particles into a container using their Monte Carlo model. A similar simulation was performed with the contact model as shown in Fig. 14. Five hundred disks were randomly positioned in a rigid container, and then allowed to fall until a stable state was obtained. A plot of the system’s potential energy as a fraction of the global minimum potential energy is given in Fig. 15. (The minimum global potential energy is the energy of the hexagonal, close-packed configuration.) We can see that after approximately 0.5 s (i.e., 500 steps using a time step of 0.001 s), the changes in energy become negligible. This is in good accordance with the results obtained by Rosato et al.

62

Time - 0.0 ms

Time - 7.0 ms t =

6.21~1s

Time - 14.0 ms

Time - 19.0 ms

Time - 26.0 ms

t = 26.2 ms Time - 33.0 ms

t =32.311x+ Time - 39.0 ms

t =38.5ms Time - 46.0 ms

t = 46.2 ms Time - 53.0 ms

t = 52.4 rns Time - 59.0 ms

Fig. 12. Snapshots of a toppling dominoes simulation, using plastic commercial dominoes with spacing A/h _ 1.7.

The area fraction occupied by the disks at the end of the simulation was computed by integrating the area under the top disks. A value of 0.80 was obtained, as opposed to 0.805 in the Monte Carlo simulation and 0.9069 for a perfect hexagonal packing. Both models give values that fall within Berryman’s bounds of 0.80 to 0.84 [28] (calculated using nearest-neighbour radii probability distribution functions). The area fraction is expected to increase slightly for larger packing heights, since the model includes gravity forces and, therefore there is greater floor stress for a high stack than foi a low stack.

4.3. Bin flows Two bin flow situations are simulated here using particles of various shapes and sizes. The first simulation shows rough grains discharging from a square bin. The second simulation uses a similar bin but discharges elongated (stick-like) grains.

t = 60.1 m.s

Fig. 13. High-speed photographs of toppling dominoes, using plastic commercial dominoes with spacing A//I = 1.7. Courtesy of W.J. Stronge.

Time = 0.0 s

Time = 0.2 s

Time = 0.6 s

Fig. 14. Three selected snapshots for the pouring of 500 equalsize disk-like frictionless particles using an average coefficient of restitution of 0.8.

63

0.6

0.4 Tie

(s)

Fig. 15. Fraction of the minimum energy (energy of hexagonal, close-packed configuration) with time for the pouring of 500 equal-size disks into a container. The dashed line represents the curve obtained by Rosato et al. as shown in ref. 3.

Time = 0.00 s

Time = 0.49 s

Time = 0.04 s

Time = 0.57 s

Time = 0.10s

Time = 0.68 s

Time = 0.17 s

Time = 0.79 s

Time - 0.29 s

Time = 0.95 s

Fig. 16. Selected snapshots from bin flow simulation using granular particles of random shapes and sizes.

4.3.1. Square bin with rough grains Figure 16 shows snapshots of the simulation using a square bin having an opening of approximately five particle diameters. The inter-particle and wall-particle coefficients of friction were given a value of 0.2. Initially, the opening of the bin is assumed to be closed. The initial configuration of the grains is obtained by first randomly positioning the grains in the bin and then allowing them to settle. The bin opening is then

assumed to open instantaneously at the start of the simulation. A study of the simulation shows the formation of temporary bridges above the bin opening. Each of these bridges collapses quickly in response to impacts with surrounding particles. In Fig. 16, the snapshot of the system at time=0.17 s shows one temporary bridge. The non-uniformity of the flow of particles through the hopper opening is indicative of its unsteady nature. These results are in good qualitative agreement with computational and experimental results obtained by Walton [29] and by Tanaka et al. [30]. Indeed, the formation and collapse of bridges in bin and hopper flows is a phenomenon widely observed in engineering practice. This phenomenon is also present in the next set of simulations. Quantitative comparison will be sought in further work. 4.3.2. Square bin with elongated grains A similar square bin is used here to discharge elongated (long, thin) grains. The initial configuration is obtained as before, assuming the bin opening is closed during the initialization. The coefficients of inter-particle and particle-wall friction are taken to be 0.4 and 0.3, respectively. The corresponding coefficients of restitution are chosen to be 0.3 and 0.5. As can be observed from Fig. 17, several temporary bridges are formed and then collapse throughout the simulation. The alignment of the grains in the bridges and through the opening should also be noted, as this represents a realistic behaviour for elongated grains (see for example ref. 31). The apparently unstable grouping of grains depends on their shape and the magnitude of the tangential friction. 4.4. Discussion on eficiency Presently, the computer implementation is not optimized. The domino-like simulations required several minutes of computation on a persona1 NeXTstation. The bin flow simulation required about 2 h of computation. It is believedthat optimization of the computer code will increase the speed of the computation by at least a factor of 3. It is difficult for the authors to make a detailed comparison of run times between the simulations using the contact model and those using a spring and dashpot model, as only a limited number of simple examples were simulated using the latter approach [2]. However, in these examples, the speed of the computation of the simulations using the contact model is about an order of magnitude faster than that for the simulations obtained using a spring and dashpot model. Once the computer implementation of the contact model is optimized, extensive large-scale simulations

nn

Time = 1.30 s

Time = 1.70 s

Time = 0.75 s

Time = 1 .OO s

fine

Time = 3.00 s

= 2.30 s

Fig. 17. Selected snapshots from square bin flow simulation using elongated grains.

will be performed, the results of which shall be reported later.

5. Conclusions

Given the need to simulate the behaviour of granular particles efficiently, this paper describes a tractable computational model that enables complex simulations of intricate particle shapes to be made. In particular, the paper ‘describes a small-scale contact model for granular particles of various shapes, sizes and material properties moving in two dimensions. Some results obtained with the model are presented and show that the contact model offers significant computational savings while apparently mimicking real-world behaviour, and being in good agreement with existing models and analytical solutions. For certain classes of applications, simulations show an order of magnitude improvement in speed over other models tested.

Two key features of the contact model are the polar representation of the particle shape, and the friction modelling. The polar shape representation greatly simplifies the contact detection process, while being sufficiently flexible to model common granular particle shapes. The friction modelling allows a high level of accuracy to be achieved, at minimal computational expense. The overall performance of the contact model varies according to the application. For example, owing to the binary impact assumption, simulations of systems containing dense, slow-moving particles are likely to predict motion slower than actually occurs. Although particle velocities are updated correctly, there may be a loss of total displacement because the algorithm leaves particles in their original positions until free to move a step without overlap occurring. Further validation is required to evaluate more precisely the effect of this binary assumption on dense systems. Further research includes the investigation of the behaviour of arbitrary granular materials in larger-scale systems, in areas such as hopper and chute flow, the discharge ofbins and silos and vibratory screening. Our numerical model can allow contact and wall stresses to be predicted, and this feature will be added to our model. Experiments will be performed to validate these simulations. This requires a sophisticated measuring system, as the collection of meaningful data for comparison purposes is a difficult task. We also plan to examine how to extend the computational methods described in this paper to threedimensional problems. Satisfactory computer modelling of granular flow must extend to three dimensions, because many flow phenomena are essentially threedimensional in character. The difficulties of such .3D modelling are non-trivial, e.g., the analysis of frictional losses is a complex task, since changes in velocity are not planar, but faster computing equipment and better algorithms may bring this goal within reach.

References T.G. Drake and R.L. Shreve, .Z. Rheol., 30 (1986) 981-993. C. Hogue, Ballistic separation, First year report of Ph.D., University of Cambridge, UK, 1990. A. Rosato, F. Prinz, K.J. Standburg and R. Swenson, Powder Technol., 49 (1986) 59-69. C. Hogue and D.E. Newland, Znt. Co@ Bulk Material Handling London, UK, Oct. 1991, pp. 217-224. S.B. Savage and C.K.K. Lun, J. Fluid Mech., 189 (1988) 311-335. J.T. Jenkins, in R.J. Kops and A.A. Lacey (eds.), Non-classical Continuum Mechanics: Abstract, Techniques and Applications, Cambridge University Press, Cambridge, UK, 1987, pp. 213-224.

65 7 N.L. Ackermann and H.H. Shen, in J.T. Jenkins and M. Satake (eds.), Mechanics of Granular Material: New Modeh and Constitutive Relations, Elsevier, Amsterdam, 1983, pp. 295-304. I3 S.B. Savage, J. Fluid Meek., 194 (1988) 457-478. 9 P.A. Cundall and O.D.L. Strack, Giotecknique, 29 (1979) 47-65. 10 3. Ghaboussi and R. Barbosa, ht. J. Numerical Anal. Methods Geomeck., 14 (1990) 451-472. 11 J.M. Ting and B.T. Corkum, fioc. 3rd ht. Con6 Computing in Civil Engineering, Vancouver, BC, Canada, Aug. 1988, 587-594. 12 O.R. Walton, in J.T. Jenkins and M. Satake (eds.), Mechanics of Granular Materials: New Models and Constitutive Relations, Elsevier, Amsterdam, 1983, pp. 327-338. 13 P.K. Haff and B.T. Werner, Computer simulation of the mechanical sorting of grains, Powder Tecknol., 48 (1986) 239-245. 14 Y. Tsuji, T. Tanaka and T. Ishida, Powder Tecknol., 71 (1992) 239. 15 G.C. Barker and A. Mehta, Pkys. Rev. A, 45 (1992) 3435-3446. 16 G.W. Hawkins, in J.T. Jenkins and M. Satake (eds.),Meckanics of Granular Materiab New Models and Constitutive Relations, Elsevier, Amsterdam, 1983, pp. 305-312. 17 P. Goldreich and S. Tremaine, Icarus, 34 (1978) 227-239. 18 B.T. Corkum and M.J. Ting, The discrete element method in geotechnical engineering, University of Toronto, Internal Rep., 1986. 19 M.P. Allen and D.J. Tildesley, Computer Simulations ofliquids, Oxford Science Publications, Oxford, UK, 1990. 20 C. Hogue, Computer modelling of the motion of granular particles, Ph.D. 7%esis, University of Cambridge, UK, 1993. 21 W.J. Stronge, Proc. R Sot. London, Ser. A, 431(1990) 169-181. 22 Y. Wang and M.T. Mason, J. Appl. Meek., 59 (1992) 635-642. 1993. 23 W.J. Stronge, personal communication, 24 J.P. Tremblay and R. Manohar, Discrete Mathematical Structures witk Applications to Computer Science, McGraw-Hill Computer Science Series, McGraw-Hill, New York, 1975. 25 S. Warshall, J. Assoc. Comput. Mach., 9 (1962) 11-12. 26 W.J. Stronge, Proc. R. Sot. London, Ser. A, 409 (1987) 199-208. 27 W.J. Stronge and D. Shu, Proc. R. Sot. London, Ser. A, 418 (1988) 155-163. 28 J. G. Berryman, in M. Shahinpoor (ed.), Advances in the Mechanics and the Flow of Granular Materials, Trans. Tech, Aedermannsdorf, Switzerland, 1983, pp. l-15. 29 O.R. Walton, 4th ht. Conf Numerical Methods in Geomeckanics, Edmonton, Alta, Canada, 1982, pp. 1261-1268. 30 T. Tanaka, Y. Kajiwara and T. Inada, Trans. Iron Steel Inst. Jpn., 28 (1988) 907-915. 31 G.W. Baxter and R.P. Behringer, Pkysica D, 5I (1991) 465-471. 32 J.T. Jenkins and M. Satake (eds.), Mechanics of Granular Materials: New Models and Constitutive Relations, Elsevier, Amsterdam, 1983.

Appendix To illustrate how the equations in section 3.6 for the impulses using Stronge’s theory are obtained, the case where sticking occurs in restitution is given here in detail. The sign convention used is such that the work done on the particle is positive during restitution and negative during compression.

Recall that B,, B, and B, depend on the inertia of the system and are given by

(A-1)

B,~+~+$+$ 2 mi

mj

Ii

4

(A-2) (A-3)

where Ii is the moment of inertia of body i at the centre of mass, which is the component of Ii Rerpendicular to the plane of motion. Let &>O be an initial sliding speed such that sliding halts during restitution. For clarity, s (s = sgn(S,) = 1 in this case) is not included in this derivation, but is included in the general equations given earlier in the paper. At the end of the compression stage, C=O, and from the integration of equation (19), the normal impulse and the work done by this component of impulse (which corresponds to the change in kinetic energy during compression) can be written as &= -C,(B,+/&)-l

(A-4)

21+‘Y=-(B2+~,)p$

(A-5)

During the restitution period, the impulse and the work done when sliding first halts are .iy=S,,(B3+p&-1

2&-

l&J =

(A-6) (B2 + /LB~)(&&)~

(A-7)

where WY represents the work done until sliding first halts. At this impulse, the normal component of the relative velocity is @Y) =

V32

+

PM&&)

(A-8)

Since sticking now occurs, the normal component the velocity at the end of restitution is

of

(&-i)y> where the coefficient of friction is now set to p =B3/ since for values of p> BJB, this velocity will not change. Thus the work done during this last period and the work done during the restitution period are

B,,

+2(Bz+~,)(~~-~,)(~y-~~)

(A-10)

66

This inequality is always valid in this case, since by substituting eqns. (A-5), (A-7) and the definition of e, (eqn. (16)) into the inequality, it reduces to the following energy constraint:

&-Iiy)’ + 2(B2 + /.&)(I$-&)(&&)

(A-11)

Using the definition of the energetic coefficient of restitution (16) and eqns. (A-11) and (A-5), the following quadratic equation is obtained for the terminal normal impulse: (&-&)”

l&G I&

(A-15)

which is always valid, since the work is negative in compression and positive in restitution. Thus the terminal normal impulse for the case where sliding halts during restitution is given by

+2(&-&4(~2+cLB3)&&4

+(B2+pB3)(Py-~y)2-e2*

(B2+j.fB3)I+=0 (A-12)

By writing this quadratic equation in the form ax2 + 2bx + c = 0, where x represents (pY-IsY), we find that a

-x=-1*

b

s

The total work done by the normal impulse during the impact can be derived from the definition of the energetic coefficient as

J

Using inequalities (24) and (27), we know that a > 0, b>O and n>O; thus (a/b)x>O. Therefore for eqn. (A12) to have only one valid solution requires that c G 0, and since (& + f13) > 0 this condition can be written as (Fy-Fy)‘~e2*i%

(A-16)

(A-14)

2FPY=2W& -e”,) = - (1 -e2,)(B2+pB3pY

(A-17) (A-18)

which shows indeed that e * = 1 represents an elastic collision. A similar analysis is used to obtain the tangential terminal impulse, and so the details are not given here.