Implementation of the discrete element method in 3D

The study of wall friction with the aid of the discrete element (DE) technique cannot be imagined without efficient and reliable computer code. The three-dimensional discrete element simulations presented in this thesis have been conducted with the aid of a computer code implemented by Goda [93] in C++. The computational algorithm implemented is based on Cundall and Strack’s [79] time stepping explicit scheme where the modeling of the motion and interaction of arbitrarily-sized spherical particles involves the execution of many thou-sands of time steps. At each step, Newton’s second law is integrated to obtain updated veloci-ties and new positions from the forces acting on the particles. The contact forces governing both the translational and rotational motion of particles are computed from normal and tan-gential contact laws. In order to develop a simulation tool which is competitive with commer-cial discrete element software the author of this thesis made his computer code suitable for input data processing, numerical model specification, graphical visualization, output data pro-cessing and graphical/text file creation. The main aim of this chapter is to give insight into the implementation process of this general purpose simulation tool.

Basic assumptions of the computational algorithm implemented are as follows. (a) Particles are rigid bodies and are permitted to overlap with one another (soft contact ap-proach). (b) Contacts occur over a vanishing small area (point contact). (c) The amount of the overlap is related to the contact force via the force-overlap (force-deformation) contact law.

(d) The time step chosen is so small that the forces acting on any particle are determined ex-clusively by its interaction with the particles and walls it is in contact with. (e) The wall ele-ments do not interact with one other, they only interact with particles. The most important characteristics of walls (wall elements) are as follows. (a) Planar walls with three or four edg-es are used (see Fig. 5.1). (b) Only one side of the wall is considered to be active. (The unit

normal vector n of the planar wall identifies the active side of the wall which can handle the ball\wall interaction.) (c) Walls may be activated and deactivated at any time instant of simu-lation. (d) The position of the walls is fixed. (Dynamic walls and walls constructed from par-ticles are not realized yet.) (e) The active side of the wall is defined by the numbering (label-ing) of vertices (right hand rule).

Fig. 5.1. Planar walls with three and four edges

The particles and walls interact with one another via the forces arisen at the contacts.

However the position of the walls is specified by the user i.e. forces acting on a wall do not influence its position. The contact force arisen between two particles or a particle and a wall is decomposed into normal component acting in the direction of the unit normal vector and shear (tangential) component acting in the common contact plane. The contact laws relate these two components of contact force to the corresponding components of the relative dis-placement via the normal and shear stiffness of the contact.

The current version of the program operates in menu-driven mode on any IBM-compatible computers running Microsoft Windows operating system. Graphics facilities im-plemented in the code allow wireframe and solid-rendered plots to be generated at any time instant (see Figs. 5.2 and 5.3). The graphical user interface allows the user to view the simula-tion space or its certain part at any stage of the simulasimula-tion. The model can be magnified on the screen for optimum viewing and all the plotted output can be saved in bitmap files. Addition-ally, simulation results may be saved in output files for later visualization or processing. The smallest possible model that may be analyzed with the program consists of a single particle.

Most problems, however, are defined by assemblies of hundreds or thousands of particles. If two particles are collided, then a point contact forms. Each particle can touch several other particles while the number of contacts can change dynamically as the particles move. The user can specify the numerical model via a data (input) file. This file is created by the user and contains commands, geometrical sizes, parameter values, etc. needed for the definition of dis-crete element model. The mapping of 3D simulation space onto 2D is realized using both par-allel and perspective projection. From particles color or grayscale-rendered while from walls solid-rendered, wireframe or combined plots can be generated at any stage of the simulation (see Fig. 5.3). Techniques and methods of computer graphics used here can be found in Szirmay-Kalos’s book [86]. Main advantage of this approach is that it yields a standalone computer code where it is not needed to use external software for graphical visualization.

Normal components of the contact forces are displayed graphically through force lines (see Figs. 5.4 and 5.5). The widths of the lines are proportional to the magnitudes of the forces. In case of particle-particle contact, the color of force line is blue and their length is equal to the sum of the radii of the contacting particles. For particle-wall contact the force lines have red color and their length is equal to the distance between the particle center and the contact point.

Fig. 5.2. Color-rendered plot of particles and walls. Snapshots from the discrete element soft-ware developed.

(a) (b) (c)

Fig. 5.3. (a) Wireframe, (b) solid-rendered and (c) combined plot of walls. Snapshots from the discrete element software developed.

(a) (b)

Fig. 5.4. Plots of contact forces (a) with and (b) without particles. Snapshots from the discrete element software developed.

(a) (b) (c)

Fig. 5.5. (a) Free falling spherical balls approaching a wall element, (b) contact forces arisen at time instant t₁, (c) steady-state contact forces. Snapshots from the discrete element

soft-ware developed.

In order to detect the contacts the simulation space is divided into 3D (cubic) cells (see Fig. 5.6). This method is known in the literature as link-cell or grid method. The size of the grid cells is equal to or slightly larger than the diameter of the largest particle. If the center point of a particle is within a given cell then the particle is assigned to the cell. This approach allows us to identify the interacting particles through the investigation of particles assigned to the neighboring cells only.

Fig. 5.6. Division of the simulation space into cubic cells. Snapshot from the discrete element software developed. One part of cells is visualized only.

To set up a model and run a simulation it is required to specify the particles (In the current version of the program, creation of new particles and/or removing of already existing particles may be realized at any time instant of the simulation.), the normal and tangential contact law, the material properties and the boundary conditions (wall elements). The

re-sponse of the assembly is studied by adjusting the parameters at particle and contact level.

The constitutive behavior is then derived automatically from the contact model and the mate-rial properties. When the particle interacts with another particle (particle-particle contact) or with wall element (particle-wall contact) contact force arises. The contact occurs at a single point through which the contact force acts. The role of the contact model is to describe the physical behavior occurring at a contact. Although a particle assembly may exhibit complex nonlinear constitutive behavior, it is achieved through the use of relatively simple contact laws incorporating stiffness properties (normal and shear stiffness). When the shear compo-nent of the contact force reaches the maximum allowable shear contact force (defined as the coefficient of friction of the contacting objects multiplied by the magnitude of the clamping force) the objects in contact slide relative to one another (sliding phase).

In order to integrate the equations of motion governing the translational and the rota-tional motion of particles the method of central differences (Feynman-Newton method) with a fixed time step t is used. Most important advantage of this method is the fastness and the great accuracy. The method can be summarized as follows: the position of a particle is com-puted at the end of a time step by using the velocity in the middle of a t interval.

t

where x(t) is the particle position at time t. In order to find the velocity halfway through the t interval, the following equation is used:

where a is the acceleration of the particle. To start the calculation, the following equation can be used:

In order to compute the contact forces simple contact laws are implemented. When two particles or a particle and a wall are in contact, the normal component of the contact force is given by the Hertz-Kuwabara-Kono equation [87] derived for viscoelastic materials (nonline-ar spring resulting in repulsive (Hertz) force connected p(nonline-arallel with viscous dissipation re-sulting in velocity-dependent damping): be-tween the contacting entities, n is a damping parameter characterizing the inelastic behavior and vn is the normal component of the relative velocity defined at the contact point. As report-ed by Schäfer et al. [87] this normal contact law gives velocity dependent collision time and decreasing coefficient of normal restitution (ratio of the relative normal post-collisional and pre-collisional (initial) velocity) with increasing initial relative normal velocity in agreement with experimental results. Since the cohesion is not taken into consideration in the current version of the computer code, the normal contact force computed is always repulsive. As the reduced mass (m_eff ) is not involved in the damping term of Eq. 5.6 the contact model is una-ble to take into consideration the mass dependence of the coefficient of normal restitution as it is concluded by Schäfer et al. [87]. The tangential contact law implemented takes into consid-eration both the tangential elasticity and the Coulomb friction law. The former allows coeffi-cient of tangential restitution (ratio of the relative tangential post-collisional and

pre-collisional velocity) to be not only positive but also negative while the latter limits the magni-tude of the tangential (friction) force. The negative coefficient of tangential restitution indi-cates that the direction of the tangential velocity is changed during the collision. The tangen-tial contact law implemented may be written as ([87])

t and v_t is the tangential component of the relative velocity defined at the contact point. At this point, however, it must be mentioned that the damping constant _t has no clear physical in-terpretation [87]. In order to obtain a non-zero friction force at rest (when vt=0), a “virtual”

tangential spring (”Cundall spring” [79]) is located between the contacting entities when two particles or a particle and a wall start to touch one another. The total elongation of this spring is set to zero at the initiation of the contact and then it is integrated over the entire collision time:

When contacting entities are no longer in contact with one other the “virtual” spring is removed. According to the Coulomb criterion, the maximum possible value of the frictional force, however, is proportional to the normal (clamping) force multiplied by the dynamic co-efficient of friction . Consequently, the tangential contact law may be written as

)

Fn/kt that is in accordance with the Coulomb criterion. As it is pointed out by Schäfer et al.

[87] the normal and tangential contact laws used here give results being in agreement with the experimental results as long as the model parameters are selected with care. [87], among oth-ers, is concerned in depth with the calibration of model parameters too. The experiments of Foerster et al. ([91]) on cellulose acetate spheres with diameter of 6 mm are used for this pur-pose. According to Schäfer et al. [87] the non-linear normal stiffness connected to the elastic properties and to the radii of the spheres through

eff connected to the radii of the spheres and the two coefficients of bulk viscosity [87]. For Foer-ster et al.’s cellulose acetate spheres the measured coefficient of normal restitution for initial relative normal velocities varying between 0.29 and 1.2 m/s is e_n 0.87 which results in damping constant of _n 190kg/m¹²s. As it is mentioned by Schäfer et al. [87], the ratio

n

t k

k / determines the results of an oblique collision where the numerical value of k_t is nearly 1000 times smaller than that of k_n. Additionally, it is also concluded that k_t has to be set such that for one particular impact velocity, the ratio of the duration of the tangential as well as the normal collision to be about one (t_t/t_n 1). In case of non-linear repulsive (Hertz) force the velocity-dependent duration of the normal collision is [87]

5 velocity. Using the “Cundall spring” the duration of the tangential collision is [87]





I  is the mass moment of inertia. For Foerster et al.’s cellulose acetate spheres m

In this chapter, the three-dimensional simulation tool of the author of this thesis is ap-plied for the simulation of wall friction in three-dimensional hoppers/silos. The results pre-sented here were reported in Goda and Ebert [93].

In order to design silos/hoppers properly and to avoid silo collapse the designer has to predict accurately, among others, the friction force and the pressure acting on the walls and the bottom of the hopper/silo for both filling and discharge. Although the hopper flow is in-tensively studied in the literature, the former three-dimensional DE simulations had been con-centrated on cylindrical/conical (axisymmetric) hoppers only. In contrast, Goda and Ebert [93] focused on rectangular hoppers and silos and not only the hopper flow and the spatial distribution of wall forces/pressures but also the effect of hopper/silo geometry on static and dynamic behavior were investigated quantitatively in 3D. Additionally, before 2005, the ma-jority of DE simulations studied hopper/silo flow in 2D and only a few three-dimensional simulations were available in this field. It must also be mentioned that the number of particles involved in the preceding analyses was, practically without exception, much smaller than used in the work of Goda and Ebert [93].

5.4.1. Literature results

Containers with a variety of shapes are widely used in industry for storage of granular materials. The wall stresses arisen at the end of the initial fill and during the subsequent dis-charge can be predicted on the basis of both analytical and numerical computations. Recently, two numerical methods are widely used to study the mechanical behavior of granular materi-als in containers: the continuum mechanics-based finite element method (FEM) and the dis-crete element method (DEM). The applicability of continuum mechanics for hopper flow, however, is strongly limited because it is unable to take into consideration the discrete nature of granular materials. On the other hand, the DEM is a powerful numerical method, in which the motion of each individual particle is determined on the basis of all the forces acting upon it (see the study of Cundall and Strack [79]). Contrary to the continuum techniques it simu-lates effects at particle level.

Many reports exist on hopper flow, but most of them are concerned with two-dimensional simulation and use cylindrical particles, as in the studies of Ristow and Herrmann [88], Masson and Martinez [89] and Sanad et al. [90], to simplify the mathematical