GRANMA manuscript No. (will be inserted by the editor)

A Fast Model for the Simulation of Non-round Particles

Alexander V. Potapov and Charles S. Campbell

Abstract This paper describes a new, computationally efficient model for the discrete element simulation of a certain class of non-round particles. The boundaries of the particles in this model are constructed from the circular segments of different radii in such a way that connections between these segments are continuous. As such, the model does not permit the simulation of arbitrarily shaped particles, but it does allow a wide enough variety of shapes to assess the effects of non-round shapes (in particular, particle interlocking) in an efficient maner. A direct test of the model’s performance demonstrates that the model is much more efficient than other models for non-round particles currently available and is less than two times slower than models for the same number of round particles.

1 Introduction In recent years, soft particle discrete element simulations have proven to be a powerful tool for investigating the behavior of granular systems. All of these simulations (see the reviews in [1,2]) follow individual particles as they move, rotate and interact with their neighbors. In a real granular system the forces generated when particles are in contact result in small changes in the particle shape, particularly around the points of contact. However, in the simulation, it is generally assumed for simplicity’s sake, that individual particles do not change shape. Instead their surfaces are allowed to overlap slightly and, in an approximation of the true elastic response, generate a repulsive response which is a function of that overlap. Virtually any imaginable dependence of the contact force on the overlap may be incorporated into the model, e.g. simple linear elastic, Hertzian or some other.

At present time, the majority of the soft particle discrete element simulations have been performed using round particles for both two-dimensional and three-dimensional systems [1,2]. The round shape is very easy to simulate in the sense that every point on its surface can be determined simply by knowing the position of the center and its diameter; in particular, there is no dependence on the orientation of the particle. As a result, it is easy to determine the overlap between two particles simply as the distance between their centers minus the sum of their radii. Also, it

Alexander V. Potapov and Charles S. Campbell

Department of Mechanical Engineering University of Southern California Los Angeles, Ca. 90089-1453, USA

reduces memory requirements as only the center position need be stored.

However the use of round particles has several limitations. Nearly all natural particles are not round, a quality which is shared even by artificial particles that are intended to be round (e.g. commercial glass beads). Furthermore, non-round particle shape can create a vast difference in the mechanical behavior because they do not roll easily. For example, bulk materials composed of round particles have angles of repose that are much smaller than natural materials (see, for example [3]). One solution to this problem is to not allow the particles to rotate at all, but instead interact only frictionally. Realistic angles of repose may be obtained using this technique, but it in no way clear that the resulting material accurately simulates all the effects of non-round particles. In particular, realistic particles may interlock with one another and produce behavior that cannot be approximated by a simple friction at the surface of round particles. Another technique is to introduce the special force term [4] to account for the Coulomb’s law behavior commonly associated with granular materials, but this moves the simulation even farther from real situation.

The most natural solution for this problem is to use particles that are not round. Several such model have been proposed. In general they can be divided into two classes. In the first, an analytic representation of the shape of a particle is used for which it is possible to simplify the contact and overlap determinations. For example, algorithms have been developed for particles with elliptic (e.g. [5–8]) or superquadric (e.g. [9,10]) shapes (the superquadric particle shape is described by the same equation as elliptic one with the only difference is that the degree is above two). Such particle can interlock like real materials, but they do not resemble natural materials. Also, aspect ratios for elliptical particles must be significantly different from unity in order to obtain realistic angles of repose, which may present other problems as the majority of natural granular materials have aspect ratios that are close to unity. The second type of model assumes that the particles are in the shapes of arbitrary convex polygons [11–15]. (Such a model forms the basis of the fracture simulations presented in [16–19].) More realistic looking particles may be created in this manner, but at tremendous computational costs. Polygonal simulations can take up to an order of magnitude more computer time than their round counterparts. For analytically represented particle shapes, (elliptic or superquadric cases) the majority of the computer time is spent in iterative solutions of the non-linear

Fig. 1. The picture of a typical oval. The boundary of this oval consists of four arches drawn around points O_{1}, O_{2}, O_{3} and O_{4}, arches drawn around points O_{1} and O_{3} have the same radius and angles α_{1} and α_{3} are the same. The same is true about the arches drawn around points O_{2} and O_{4} and angles α_{2} and α_{4}. Points O_{1} and O_{3} are situated in such a way that line connecting these points pass through the center of an oval, the same is true about the line connecting the points O_{2} and O_{4}.

equations used to calculate the particle overlap. In the polygonal simulations, the majority of the computer time is spent calculating intersections of the sides of contacting polygons which is a necessary part of the overlap determination.

2 The Description of the Model

Here we propose a soft particle discrete element model for the simulation of non-round particles which exhibits a performance that is comparable to the performance of round particle models and is able to accommodate a much wider variety of shapes (and more realistic shapes) than the elliptic or superquadric models. The idea for this model was born from the way the ellipses are sometimes approximated in technical drawings; when appropriate elliptical templates are not available, an oval with the same aspect ratio is often used to approximate the ellipse. An oval is geometric figure whose boundary is determined by four circular arches of two different radii which are joined together in a continuous way (i.e. the first derivative is continuous between the arches) and thus can be drawn with compasses and other standard tools used to draw circles. However, it is possible to generate more complex shapes, many nearly polygonal in shape, from circular arches. At the same time, the determination of contact and overlap between two circular arches segments is nearly identical to the same procedures for circles. Thus a simulation of particles generated from circular arches will allow a variety of

Fig. 2. This picture illustrates the procedure of the calculation of the contact parameters between a pair of ovals. After the contacting arches are determined (arches corresponding to angles α_{1} and α_{2}), the calculation of the contact parameters can be carried out in the same way as it is done for the round particles.

shapes to be studied at only a small computational cost beyond that of round particles.

To illustrate this technique, consider the typical oval presented in Fig. 1. The boundary consists of four arches drawn around points O_{1}, O_{2}, O_{3} and O_{4}. In a classical oval, the arches drawn around points O_{1} and O_{3} have the same radius (we shall call this radius “large”) and the angles α_{1} and α_{3} are identical. The same is true about the arches drawn around points O_{2} and O_{4} (we shall call this radius “small”) and the angles α_{2} and α_{4}. Points O_{1} and O_{3} are situated in such a way that line connecting these points passes through the center of an oval, as does the line connecting the points O_{2} and O_{4}.

Now suppose that we have a granular medium consisting of oval shaped particles which will be simulated using an otherwise standard soft-particle discrete element model. To do this, we must at every time step find the particles which are overlapping, calculate the distance of this overlap and determine the center of this overlap in order to know the point at which contact forces are to be applied to the contacting particles. Let us suppose that we have a contact between two ovals as is shown in Fig. 2. Any possible contact must occur between the arch of one oval and the arch of the other oval. This type of contact can be treated in the same way as the contact of the two circular particles. Thus the overlap distance is simply the sum of the two radii of the contacting arches minus distance between the points around which these arches are drawn (i.e. distance between the points O_{1} and O_{2} in Figure 2). The center of the overlap lies in the line connecting O_{1} and O_{2} midway between the points of intersection of this line and the two arches (point C in the Figure 2).

Compared to the case of round particles, the only additional operation which is required here is to determine

Fig. 3. Examples of the construction of the particles out of the arches continuously joined to each other (the first derivative is the same on the both sides of the arch to arch transition point):

(a) approximation of a triangle, (b) approximation of a square,

(c) approximation of a pentagon and (d) approximation of a heptagon.

which of the arches may be involved in the contact. This can be carried out in very simple and effective way. Note that two ovals can have the contact along a pair of arches if and only if the line connecting the points around which these arches are drawn lies within the limits of the angle of both arches (in Figure 2 this requires determining that the line O_{1}O_{2} lies between angles α_{1} and α_{2}). This is easy to check by calculating the dot products of the unit vector connecting the points around which the arches are drawn with unit vectors of the bisectors of the angles of the arches. Plus, the dot product of the vector connecting the centers of gravity of the ovals and the vector connecting the points around which the arches are drawn should be calculated to ensure, by its sign, that the potential contacting arches are oriented towards rather than away from each other. The calculation of these three dot products is an additional computational step with respect to the case of the round particles. However, it is not necessary to perform all three calculations at each time step. Since the particles do not drastically change their orientation during a single time step, in the vast majority of cases, it will only be necessary to check at each time step whether the same two arches remain in contact. This means the calculation of only the first two dot products, which is rather fast and straightforward procedure. Also, it is quite likely that any two arches that are in contact, were in contact during the previous time step. Consequently, much computer time can be saved by simply verifying that the contact continues on subsequent time steps. The calculation of the contacts of the ovals with the straight boundary walls is performed in much the same way.

The model described above is easily extended beyond the case of ovals. Ovals, like ellipses, have the disadvantage of requiring relatively large aspect ratios to obtain realistic values of the angles of repose. Much more complex particles may be constructed in the same way the ovals are constructed. If the boundary of a particle consists of any number of arches which are continuously connected together, one can see that the algorithm for determining contact as described above still holds true. Some examples of this sort of particles are presented in Fig. 3. For these particles, a large radius is used to construct the “sides” of the particles and a smaller radius is used to connected them smoothly together at the “vertices.” Clearly, one can easily construct an approximation to a triangle, square, pentagon and so on, and one can get very close to an actual polygonal shape as one wishes by making the corresponding radius of curvature of the arches that form the sides very large and those that form the corners, very small. (This process is limited only by roundoff errors in the determination of the overlap. I.e. as the overlap is computed by taking differences, high accuracy requires that its magnitude be large enough relative to the radius to be represented by several number significant digits in its floating point representation within the computer.) If these particles are constructed of the arches with only two values of radii as it was done for the oval and all of the shapes in Fig. 3, no additional memory (beyond that required for an oval) need be used for the storage of the information about position and orientation of the particles.

One additional advantage and one additional disadvantage of the model proposed here should be mentioned. The advantage is that this model has a fixed value of the radius of curvature at every point of the surface of the particles that is bounded by continuously joined arches. This makes it possible to use the Hertzian contact model (which describes the contact between particles with arbitrary, but finite and non-zero, radii of curvature) without any additional assumptions - which would be a questionable for the case of the angular particles. The disadvantage of the model is that for the moment we can not see its simple three-dimensional implementation.

3 Demonstration of the Performance of the Model It is clear from the material presented above that the performance of the proposed model should be only marginally worse than a round particle model as only a few additional calculations are required. A direct test of this fact is described in the following.

For the test, we have chosen the problem of settling the particles into a rectangular box under the influence of the gravity. In the first case, the particles are round, in the second case, the particles are the same number of quasitriangles (Figure 3a), in the third case the particles are quasi-squares (Figure 3b), in the forth case, particles are of superquadric shape with superquadric degree 3.0 which are simulated by the technique described in [10] and finally, in the fifth case particles are square and simulated by the algorithm described in [12]. Linear spring contact laws are assumed for all cases except for the case of square

Fig. 4. Results of the simulation of gravity settling round particles. This simulation has been performed in order to compare the efficiencies of the models. (a) Initial positions of the particles; (b) positions of the particles at the end of the simulation.

Table 1. CPU time spent on SPARC-20 for performimg of 4.69 · 10^{4} time steps

Particle shape CPU time (seconds)
Round  29
Quasi-triangle 48
Quasi-square 51
Superquadric 347
Square  114

particles for which the normal contact force was based on the area of overlap (see [12]). The geometry of the boundary walls, the particle masses, particle densities, the time step and the coefficient of the restitution of the particles were the same in all cases. The models were run on the same Sparcstation 20 computer (75 MHZ processor), for the same period of simulation time (4.69·10^{4} time steps), and the amount of CPU time expended was determined and compared. The same contact search routine for near neighbors [12] was used in all simulations.

The initial positions of the particles are presented in Figure 4a for the case of round particles, in Figure 5a for the case of quasi-triangular particles, in Figure 6a for quasi-square particles, in Figure 7a for the superquadric particles and in Figure 8a for square particles. Figures 4b, 5b and 6b, 7b and 8b depict the final positions of the particles for these five cases.

The results of the simulation of these cases on a SPARC-20 workstation are presented in Table 1. These results show the CPU difference of about 66% between the case of round particles and case of the quasi-triangular partic-

Fig. 5. Results of the simulation of the gravity settling of the quasi-triangular particles in a rectangular box. It was performed for the same particle mass, same time stepand same time elapsed in the system (same number of time steps) as the simulation for round particles presented on Figure 4. (a) Initial positions of the particles; (b) positions of the particles at the end of the simulation.

Fig. 6. Same as Figure 5, but for quasi-square particles

les composed of arches. The case of quasi-square particles is just 4% slower than the case of quasi-triangular particles, this time difference is due to slightly longer procedure of the determination of contacting arches for quasi-square particles than for quasi-triangular particles due to larger number of arches. The case of superquadric particles was more than one order of magnitude slower than the case of

Fig. 7. Results of the simulation of gravity settling of the same number of superquadric particles [9] as on Figure 4. It was performed for the same particle mass, same time step and same time elapsed in the system (same number of time steps) as the simulation for round particles presented on Figure 4. (a) Initial positions of the particles; (b) positions of the particles at the end of the simulation

Fig. 8. Results of the simulation of gravity settling of the same number of square particles as on Figure 4. It was performed for the same particle mass, same time step and same time elapsed in the system (same number of time steps) as the simulation for non-round particles presented on Figure 4. (a) Initial positions of the particles; (b) positions of the particles at the end of the simulation.

round particles. (It is hard to do exact performance comparisons for this case since the procedure is iterational and “exact” solution is never reached. We used four iterations per time step here, different number of iterations may produce different results.) Finally, the model with square particles was almost four times slower than the model for round particles and thus more than two times slower that the model proposed in this paper. Thus, the model presented here easily outperforms every other model for non-round particles available at the present time.

4 Conclusions In this paper, a model has been proposed for a soft-particle discrete element simulation of the motion of non-round particles. This model is based on the construction of a particle boundary out of a set of arches which are connected to each other so that the first derivative at the point of transition from one arch to another is the same from both sides. The model allows us to vary considerably the degree of roughness of a particle surface, and it performs approximately 60% slower than a model for the round particles, thus easily outperforming any other model currently available for non-round particles.

The model is capable of efficiently simulating nearly polygonal particles. The main use of this type of simulation is not necessarily to simulate the actual shapes of the particles in naturally occurring materials, but to provide a way of adding the effects of non-round shapes (in particular, particle interlocking) into existing simulations.

References

1. Campbell, C.S., Computer simulation of rapid granular flows, Proc. 10^{th} US National Congress of Applied Mechanics, Austin Texas, June 1986, ASME, New York, p.327-338.

2. Campbell, C.S., Computer simulation of powder flows, to appear in Powder Techology Handbook, Second Edition. (Gotoh et al. eds) Dekker, New York, 1997, 777-793.

3. Poschel, T. & Buchholtz, V., Static Friction Phenomena in Granular Materials: Coulomb Law versus Particle Geometry, Physical Review Letters 71 (1993), p. 3963-3966.

4. Cundall, P.A. & Strack, O.D.L., A discrete numerical model for granular assemblies, Geotechnique 29 (1979), p.47-65.

5. Rothenberg, L. & Bathurst, R.J., Numerical Simulation of Idealized Granular Assemblies with Plane Elliptical Particles, Computers and Geotechnics 11 (1991), p.315-329.

6. Rothenburg, L. & Bathurst, R.J., Micromechanical features of granular assemblies with plane elliptical particles, Geotechnique 1 (1992), p.79-95.

7. Ting, J.M., A robust algorithm for ellipse-based discrete element modeling of granular materials, Computers and Geotechnics 13 (1992), p.175-186.

8. Kohring, G.A., Computer simulations of sintering via granular dynamics, Physica A 195 (1993), p.1-11.

9. Williams, J.R. & Pentland, A.P., Superquadric and Modal Dynamics for Discrete Elements in Concurrent Design, Proc. 1^{st} US Conference on Discrete Element Methods, Golden Colorado, October 1989, (Eds. Mustoe, G.G.W., Henriksen, M. & Huttelmaier, H.P.)

10. Mustoe, G.G.W. & DePooter, G., A numerical model for the mechanical behavior of particulate media containing non-circular shaped particles, Powders and Grains 93, 1993, (Ed. Thornton, C.) Balkena, Rotterdam, p.421-427.

11. Hogue, C. & Newland, D.E., Efficient computer modeling of the motion of arbitrary grains, Powders and Grains 93, 1993, (Ed. Thornton, C.) Balkena, Rotterdam, p.413-419.

12. Hopkins, M.A., The Numerical Simulation of Systems of Multitudinous Polygonal Blocks, US Army Cold Regions Research and Engineering Laboratory, USACRREL Report CR 99-22, 1992.

13. Kohring, G.A., Melin, S., Puhl, H., Tillemans, H.J., Vermohlen, W., Computer simulations of critical, nonstationary granular flow through a hopper, Computer Methods in Applied Mechanics and Engineering 124 (1995), p.273-281.

14. Tillemans, H.J. & Herrmann, H.J., Simulating deformations of granular solids under shear, Physica A 217 (1995), p.261-288.

15. Poschel, T. & Buchholtz, V., Molecular Dynamics of Arbitrarily Shaped Granular Particles, Journal de Physique I 5 (1995), p.1431-1455.

16. Potapov, A.V., Hopkins, M.A. & Campbell, C.S., A Twodimensional Dynamic Simulation of Solid Fracture. Part I: Description of the Model, International Journal of Modern Physics C 6 (1995), p.371-398.

17. Potapov, A.V., Campbell, C.S. & Hopkins, M.A., A Twodimensional Dynamic Simulation of Solid Fracture. Part II: Examples, International Journal of Modern Physics C 6 (1995), p.399-425.

18. Potapov, A.V. & Campbell, C.S., A Hybrid Finite-element simulation of Solid Fracture, International Journal of Modern Physics C 7 (1996), p.155-180.

19. Potapov, A.V. & Campbell, C.S., A Three-dimensional Dynamic Simulation of Solid Fracture, International Journal of Modern Physics C 7 (1996), p.717-730.