Discrete Element Methods: Numerical Modeling of Discontinua

The case is made that continuum methods (such as finite elements) for rock and soil may be completely replaced by particle models (such as DEM) in 10–20 years. Reproducing complex behavior with continuum methods requires complex constitutive models, containing sometimes dozens of parameters and/or internal variables, such as yield-surface descriptors. The DEM model needs very few material properties because the behavioral complexity arises as an emergent property of the system. There is a lot of current interest about emergent systems in which collections of simple things behave—collectively—in complicated, and completely different, ways. For example, in rock and soil, the following behaviors appear "automatically" from a particle model of the material: continuously nonlinear stress/strain response, dilation related to mean stress, transition from brittle to ductile behavior, hysteresis and memory (nested yield surfaces), nonlinear strength envelopes, acoustic emissions, and many more. It is very complicated to reproduce all these things with a continuum constitutive law, and such laws tend to contain many ad hoc parameters and equations. However, all the phenomena noted above appear as emergent properties of a macroscopic system composed of particles—either unbonded for soil, or bonded for rock. Further, a DEM model naturally exhibits localization (fractures in a brittle solid, or shear bands in a granular material). It is difficult to capture such localization in a continuum model that uses a mesh. It is already accepted that the DEM is well suited to simulating systems such as ball mills, blocky rock and assemblies of discrete fragments. Further, the method has been used to study the micromechanics of materials such as soil at the element level. However, we suggest that in a few years the DEM should be applied in the traditional domain of continuum methods, representing complete rock & soil systems such as dams, tunnels and foundations. Clearly, this can only be done when computers become powerful enough, but it is argued that this will occur in 10–20 years. The talk will explore the reasons why all the complex behaviors of soil and rock (mentioned above) emerge at the system level, given very simple laws (and few parameters) at the micro-level. It will also be argued that the methodology enables the engineer to obtain good insight into nonlinear mechanisms that occur in deforming geomaterials. For example, a particle model of brittle rock reproduces the results expected from fracture mechanics, but the model leads to a new (and more intuitive) explanation for the behavior. Although particle models of complete geo-engineering systems must wait until computer power becomes sufficient, we should start experimenting with the methodology now. There are a number of problems to solve, such as the development of an efficient calibration procedure, before the methodology can be used routinely.

Formulations of problems of continua are based on what are termed constitutive laws. From the constitutive law and preservation laws (preservation of energy, momentum, mass, etc.), a set of governing equations defining the problem is derived. A continuum model is only an approximation of the real physical phenomena and is as good as are underlying assumptions on which it is based. The most important assumption necessary for the validity of the continuum model is the assumption that the scale of the problem or the scale of the part of the problem of practical significance is much larger than some characteristic length defined by the microstructure of the material. For elastic analysis of solids this characteristic length is one order of magnitude greater than the size of molecules making this material. Thus most of the practical problems are well represented by the continuum model. Center now on a very different set of problems such as filling a basket with pebbles. The question is how many of these pebbles can fit into a single basket? It is evident that the scale of the problem is similar to the characteristic length of the microstructure of the matter. Problems like this are of practical importance in many scientific, engineering and industrial applications. The common feature of all of them is that the representative volume of the matter is either much larger or of the similar order of magnitude as the physical problem analysed. Thus the necessary conditions for continuum based governing set of equations are not satisfied. The researchers and engineers have recognised this fact for decades. Thus the solution was sought through either experimental investigations or phenomenological analytical approximations.

The classical discrete element method (DEM) proposed by Cundall and Strack has been proved to be a versatile numerical tool, particularly suitable for the simulation of granular and particulate systems typically encountered in process engineering and geomechanics. Such systems are of an inherent discrete nature and often consist of an excessively large number of individual particles. The overall behaviour of the systems is determined by the motion of these individual particles involving interaction mainly through adhesive/frictional contact. In principle, the DEM is based on the concept that individual material elements, each usually assumed to be rigid, are considered to be separate and are (possibly) connected only along their boundaries by appropriate physically based interaction laws. Many other industrial and scientific problems are, on the other hand, characterised by a transformation from a continuum to a discontinuum state. For example, material separation and progressive failure phenomena can be found in applications such as concrete structural failure, food technology processes, rock blasting operations and fracture of ceramic or other quasi-brittle materials under high velocity impact. The problems are initially represented by a small number of discrete regions (often a single entity) prior to the deformation process. During the loading phase, the bodies are progressively damaged and modelling of the subsequent fragmentation may result in possibly three to four orders of magnitude more bodies by the end of the simulation. The overall system response is governed firstly by appropriate constitutive mechanisms which control the material separation process followed by description of the inter-element interaction forces which can be short-ranged, such as mechanical contact, and/or medium-ranged, such as attraction forces in liquid bridges, which control the subsequent motion of particles. There is currently intense research interest in the solution of this class of problems. Strategies range from continuum based approaches, such as decohesive zone models, to discontinuum driven formulations, such as discrete discontinuous analysis (DDA) methods. For problems in which interest is restricted to relatively small deformations, the use of continuum based methods may be suitable, but for situations involving large topological changes, such as the modelling of particle flow behaviour post-fracture, there are compelling advantages in employing combined finite/discrete element solution strategies. It may be argued that by modelling the continuous to discrete transformation explicitly, a physically more realistic representation is obtained, which provides advantages in allowing a constitutive description of the material that requires a minimum number of parameters, which may be identified from standard experimental tests. Over the last decade or so, such a combined solution approach has emerged as a powerful numerical tool for the simulation of large scale industrial problems with discrete/discontinuous features. The paper considers the following essential issues necessary for an effective computational implementation of a continuum/discontinuum formulation: (1) fracture criteria and propagation mechanisms within both the finite and discrete elements, together with mesh adaptivity procedures for discretisation and introduction of fracture systems; (2) detection procedures for monitoring contact between large numbers of discrete elements; (3) interaction laws governing the response of contact pairs; and (4) parallel implementation.

Particle science and technology is a rapidly developing interdisciplinary research area with its core being the understanding of the relationships between micro- and macroscopic properties of particulate/granular matter which is widely encountered but poorly understood. The macroscopic behaviour of particulate matter is controlled by the interactions between individual particles as well as interactions with surrounding gas or liquid. Understanding the microscopic mechanisms, i.e., the particle-particle and particle-fluid interactions, is key to leading to truly interdisciplinary research into particulate matter, in which scientists and engineers correlate their findings and ensure that microscopic predictions from one discipline match macroscopic results from another. This aim can be effectively achieved via particle scale research based on detailed microdynamic information such as the forces acting on and trajectories of individual particles in a considered system. In recent years, such research has been rapidly developed worldwide, mainly as a result of the rapid development of discrete particle simulation technique and computer technology. This paper will present a brief review of the work in this direction in my laboratory, "Particulate and Multiphase Processing" at the University of New South Wales (UNSW), in the past years. It covers, through demonstrative examples, the development of simulation techniques and their application to the study of particle packing and flow, transport properties and constitutive relationships of typical static or dynamic particulate systems. It is concluded that DEM is an effective method for particle scale research of particulate matter.

This paper describes application of the Finite Wall method to the problem of metering seeds singly from a bulk. The purpose of the meter is to pick-up single seeds from a reservoir and deliver them to a planting mechanism at regular intervals. The device modeled comprised a slotted bowl rotating about a horizontal axis, and enclosed within a casing with an opening through which seeds are gravity-fed from an external hopper. Seeds are captured from the bulk by semi-circular cells distributed around the bowl rim. The system represents a major challenge for DE simulation as it comprises both static and moving boundary surfaces of irregular shape interacting with non-spherical particles.

The behaviour of each object in a discrete element system is largely determined by the interaction laws specified in the context of interelement contact. The contact interaction laws are normally developed within the framework of penalty methods. The penalty method only approximately satisfies the contact displacement constraints, allowing a small amount of overlap to occur. In this context, the interaction laws describe the relationship between the overlap and the corresponding contact force between a contact pair. Normal contact laws are fairly well established for bodies with smooth boundaries, such as disks, spheres, ellipses and ellipsoids, but it is not the case for discrete bodies with non-smooth boundaries. One notable example is that there is still a lack of a simple and robust algorithm for modelling corner-comer contact of two polygonal objects. The difficulty mainly lies in the fact that it is not always clear how to determine a proper direction along which the normal contact forces should be applied due to the directional discontinuity at the corners. Similar situations also appear in contact between finite elements and frequently cause convergence problems in the nonlinear iteration procedure. In the combined finite/discrete element simulation of multi-fracturing phenomena in particular, inappropriate treatment of comer-corner contact may introduce substantial energy into the system and thus makes the computed results, unreliable. At present, most of the existing comer models employed in the discrete element method are of a heuristic nature and generally very complex, especially when extended to 3D cases. In addition, energy conservation cannot always be guaranteed in the models. In the current work, an attempt is made towards the development of a non-ambiguous corner-corner contact algorithm within an energy based theoretical framework, in which the direction, magnitude and reference position of the normal contact force are uniquely defined, thereby circumventing the difficulties mentioned above. Furthermore, it also provides a sound theoretical foundation to possibly unifying existing contact models for all types of convex objects within the framework.

Discrete element modeling has generally used simple circular shapes such as two-dimensional disks and three-dimensional spheres and not so simple polygonal and polyhedral shapes. Contact detection ranges from trivial for disks and spheres to non-trivial for polygons and polyhedra. To deal with complex shapes while circumventing the complexities of contact detection some researchers have constructed complex particles from assemblies of spheres. This strategy results in a trade-off between contact simplicity and a large increase in the number of particles, contact searches, and particle contacts.

The paper describes the torus as a primitive shape for constructing realistic particles with the Discrete Element Method (DEM). The current work is intended to complement a paper, currently in review, in which the combination of a torus and two spherical caps is shown to approximate any axisymmetric ellipsoid (Kuhn, 2002b). The range of shapes that can be constructed from torus and sphere pieces is illustrated in Figure 1. Each of these shapes was created by joining two spherical caps with a central torus. Even more shapes are possible if cylinder and planar surfaces are also included as primitives. When only spheres and tori are used, however, the shapes will be smooth and have non-zero Gaussian curvature (except, perhaps, at patches of zero measure). These features allow one to avoid the problem of analyzing contact forces at sharp edges or corners. The greatest advantage of using torus primitives is rapid contact inquiry. The author has described stable and efficient algorithms for contact inquiry between torus-sphere and torus-torus pairs. For these two cases, contact inquiry is only about three times slower than for sphere-sphere pairs, and the algorithms were found to be stable through over 500 million random trials. In the current paper, we extend previous work by presenting geometric characteristics of a general torus section. We then present an iterative algorithm for contact inquiry between a torus and a cylinder.

A new boundary condition was developed for discrete element method (DEM) simulations. The boundary is assumed to be a body of pressurized hydrostatic fluid. When an ellipsoidal particle touches the boundary plane, the interaction between the particle and the boundary is calculated based on the area and the center of the cross-section of the particle by the plane, which were determined analytically. This new boundary condition was implemented into an existing ellipsoidal discrete element code. The new code was then used to generate a specimen. The specimen contains 1170 mono-sized ellipsoids with an aspect ratio of 1.5. The specimen is isotropically compressed to 100 kPa. A triaxial compression drained test was conducted. The result showed an increase in macroscopic friction angle comparing with the values commonly obtained using periodic boundaries. Different boundary conditions have been used in the 3-D DEM simulations. The most commonly used are periodic boundaries and rigid boundaries. Uniform strain field is assumed in the periodic boundaries. Boundary movements are imposed in rigid boundaries. There are few simulations that used flexible boundaries. Additional boundary elements are used in flexible boundaries that increase the computation effort. The effect is evidential when small sample sizes are used. The macroscopic friction angle obtained from the simulations with periodic and rigid boundaries is lower than the results from triaxial tests on sands. The author expects that this low simulated macroscopic friction angle is a function of the boundaries. Boundary movements are applied in triaxial tests. Therefore, the assumed uniform strain field may not occur in triaxial tests. The author has used rigid boundaries (having the same material properties as the particles) in the DEM simulations. The friction angle did increase when the sample size is small. However, when the sample size is larger than 1000 particles, the obtained macroscopic friction angle is still lower than the experimental values. This paper will present a newly developed boundary condition. The new boundary is implemented in the DEM. The calculation of the contact state and a numerical result of a triaxial compressive drained test will be presented.

This paper presents an algorithm that determines the spatial relationships between bodies and can be used for contact detection. The process of contact detection can be divided into two phases, called neighbor search and contact resolution. Neighbor search identifies and creates lists of object "near" to the target object, usually using some approximate geometry for the objects. The geometric resolution then compares the detailed geometric representation of each object with the target object to resolve contact. Based on the details of the contact, forces or some other constraint enforcement technique is used to penalize and remove any overlap that has occurred due to the discrete time stepping algorithm. In simulations the spatial algorithm usually becomes the computational bottleneck as the number of objects increases. It is therefore critical to develop algorithms that scale well, preferably linearly, with respect to the number of bodies and polygons. In addition the algorithms should be general enough to be applicable across a wide spectrum of problems.

The discrete element method (DEM) is ideal for modeling many problems not accessible to traditional continuum-based methods such as finite difference and finite elements. DEM has the advantage of inherently capturing large non-affine deformations and the fluid—solid phase changes so commonly found in granular media. The DEM appears to derive much of its power from the kinematic freedom for the particles such that seemingly simplistic micro-scale models can replicate realistic macro-scale behavior. With the advent of large-scale computing, by which simulations of several million particles are possible, the DEM will become an important tool for geotechnical and industrial applications within the next 20 years. The advances in DEM have and will continue to mirror advances in computer hardware. Several advances in DEM technology will accompany the improvements in computer technology. The constitutive response of the DEM medium depends on micro-scale laws that define the contact mechanisms, particle size and shape, and particle distribution. The success of the method requires a correlation between the phenomenology of the DEM medium and the micro-scale laws. Most of these advances will be accomplished through numerical experiments, which at present are a popular tool for understanding granular media as a continuum. Is it reasonable to foresee a day when a continuum view is unnecessary? In this paper, this question is discussed from the standpoint of scale effects that infect all numerical approximations, whether the governing equations are derived from a continuum relationship or begin from a discrete model. The key point is that unless computing power reaches a point where there is a one-to-one correspondence between simulated particle size and actual particle size, some understanding of the DEM medium as a continuum will be needed.

Spherical discrete elements with modern contact detection algorithms are the most computationally efficient particles to use in computer simulations. However, spherical elements exhibit some natural behavior that is inconsistent with rock behavior, particularly during muck pile formation. This unnatural behavior is related to a limitation of the inter-particle friction that can be generated by the inherently smooth spherical elements. An algorithm has been developed that allows 2D and 3D explicit time integration spherical element codes (in this case DMCBLAST and DMCBLAST_3D) to determine whether a particle is in the muck pile being formed subsequent to explosive loading of the discrete elements during simulation of a rock blast. Once the particle has arrived in the muck pile, its landing/stopping behavior can be modified to better emulate the friction between rough and jagged rocks that constitute real muck-piles. Excellent correlation has been obtained between discrete element predicted muck-pile shapes and those produced by actual rock blasting and surveyed in the field. This paper presents the algorithm developed as well as 2D and 3D predictions and correlation with field data.

A contact resolution algorithm for generic shapes based on level surfaces in DEM is presented. This algorithm is only able to detect whether or not two general shapes are in contact. In the case of small particles with electrostatic forces, contact detection is not sufficient, as inter-particle forces are active even if they are not in contact. It is therefore extended here to obtain the distance between the particles, or the virtual overlap, and the initial contact point. From a user point of view, these algorithms for general irregular shapes are not very attractive. They require the solution of a system of nonlinear equations. A Newton-Raphson approach to linearise these equations leads to a series of linear equations containing the gradients and the Hessians of the level surfaces. These gradients and Hessians generally have to be obtained explicitly, either by manually deriving the function, or by some symbolic algebra package. In addition, they have to be coded and debugged by the user. These steps make the implementation of new shapes, although trivial, error-prone and time-consuming. Derivation, however, is as mechanistic process, and there is no reason why a software component being able to execute the evaluation of a function, specified by the user in the form of an algorithm, would not also be capable of providing partial derivatives of any order. Here a technique, called Automatic Differentiation (AD) is used to accomplish this. By the use of AD, the user only has to provide the code for the evaluation of the level surface. Under the hood, the software provides first and second order partial derivatives automatically to construct the linear systems for contact resolution. The use of DEM as a numerical technique for investigating granular flow problems and eventually for engineering bulk handling systems is very promising. Yet there are many scientific and numerical issues that need to be solved before these promises can become a reality. One of these problems is dealing with irregular shapes. Several approaches have been taken to treat irregular shapes. Some of them are based on polygonal representations. Polygonal objects lack curvature and consequently it is difficult to define realistic contact force models (unless a very high number of polygons is used). Overlapping spheres have also been used. Spheres do not lack curvature, but generally they provide the wrong curvature at a point and hence the wrong force. A more generic approach was taken by Cleary et al. They based their approach on level surfaces but concentrated, on hyper-ellipsoids. They provided an efficient algorithm for contact detection in the form of a quantity which depends on position and orientation of the respective particles which changes sign when in contact. For electrostatically charged particles this is not sufficient as there are inter-particle forces acting even when contact is not established and the closest distance between particles is needed and the algorithm needs modification. This problem will be tackled in section 1. A further problem of a more practical nature is that the algorithm relies on the Newton-Raphson scheme for the solution of a nonlinear system and hence the user is required to specify not only the function describing the level surface of the particles, but also the gradient and the Hessian with respect to position. The latter is necessary because the gradient already appears in the nonlinear system. Although this task is trivial in principle, it is error-prone and time-consuming. This problem will be dealt with in section 2.

Granular flow, particle mixing and powder compaction are very important processes in many industrial applications such as pharmaceutical tabletting, minerals selection, casting, packing and transport of granular materials. Flow and deformation of particles has been typically treated as a continuum process, and thus the Finite Element Method (FEM) has been widely used. In reality, however, granular powder consists of millions of tiny particles of different sizes and shapes, the conventional continuum method cannot model these discrete features. Naturally, a discrete element method can provide more realistic results using uniform spherical particles. Furthermore, powder flow and deformation modelling involve rheology, particle size distribution, fracture, and mixture of different materials. The modelling of such processes is very difficulty due to their nonlinearity. Thus, the proper modelling of powder flow and compaction is of both industrial importance and scientific interest. The main issue in discrete element modelling is the contact detection of interacting particles during the flow and deformation processes. In fact, nearly two-thirds of the computing time is used for the contact detection purpose. Thus, the efficiency of the contact detection algorithms used in the discrete element analysis is crucially important to the whole modelling procedure. A number of good contact detection algorithms have been developed such as the direct check algorithm, binary tree algorithm, contact list and grid searching algorithms. However, most of these algorithms work well for rigid spherical or triangular particles. They generally do not apply well to the non-spherical particles, or at least there is no guarantee of their efficiency. So many existing studies have focused on the flow and transport of rigid spherical particles and soft particles. Recently, there have been new developments concerning flow of non-spherical particles such as elliptical shapes and/or triangular particles. In this paper, we present a two-stage contact detection algorithm for non-spherical particles using the contact list algorithm and direct polygon contact detection. This method combines the advantages of several existing contact detection algorithms so as to find an efficient way to deal with the interactions of non-spherical particles, and simply treating the spherical particle as a degenerated case of non-spherical shapes. In this way, we can essentially deal with any particle shape. In the next sections, we will give a brief description of the two-stage contact detection algorithm, and then provide some simulation results and compare the efficiency with other contact detection procedures.

Practical geotechnical problems generally are three-dimensional in nature. The geometry of the solution or the construction works, the properties of the various rock materials and the discontinuities all possess strong three-dimensional characteristics. At present, Discrete element method (DEM) and Discontinuous deformation analysis (DDA) are mainly limited to two-dimensional analysis and there are not many three-dimensional studies. Three-dimensional analysis is however important in many real engineering problems and should be considered. In DDA and DEM analyses, the system is composed of discrete blocks. Detecting and representing contacts between blocks becomes one of the key problems of these two methods. First of all, the numerical analysis depends heavily on the contacts between blocks. Wrong detection or description of contacts may lead to wrong analysis, and breakdown of the computation may even happen. Secondly, because both DDA and DEM analyses adopt an iterative scheme to account for the large displacement and deformation, the coordinates of the blocks and contacts between blocks can change at each iterative step. The efficiency of the algorithm for detecting and representing the contacts will directly and strongly influence the efficiency of the numerical analysis of the blocks system. In fact, there are great differences in the computational time required between the use of different methods to detect and represent contacts. Detection of contacts under two-dimensional condition is a simple problem and many famous algorithms and codes are available which will not be discussed here. Detection of contacts between general 3D blocks is however much more complicated. Currently, direct determination and common plane (c-p) method are the two most popular approaches for contacts detection and representation. Direct computation is the simplest but most time-consuming approach. Compared with direct test, c-p method can greatly reduce the amount of calculation for contact detection and representation and is adopted in most of the 3D DEM and DDA programs. The computer time consumed for contact detection using c-p method is however a significant portion of the entire time for analysis. In general, there are many possible forms of 3-D blocks that make the contact detection and representation a very difficult and time-consuming task. It is however the key problem in the development of 3-D DEM and DDA programs. In this paper, a new approach for detecting and representing contacts between 3-D blocks is developed and the advantages of the new method are discussed.

Both particle packing and particle contact parameters play important roles in determining the mechanical behavior of a discrete particle assembly. One common approach to calibrate discrete element models (DEM) of rock is to hold particle packing constant while contact parameters are changed to match observed macroscopic mechanical behavior. This often leads to non-unique material calibrations, hi other words, the mechanical properties of the DEM assembly depend on the packing method. Another problem with this technique is that the parameters that describe contact interactions of particles are not directly observable in the laboratory. An opposite approach is to suggest that since we can actually observe rock texture and packing in rocks, why not mimic this in the model and keep contact parameters constant? The best method probably lies somewhere between these two end-members as suggested by mineralogical observations (i.e., not all grains have the same mechanical properties) and textural observations (i.e., our inability to capture minute details of grain shapes) of rocks. In this context, the purpose of this paper is to gain an understanding of how different particle packing and textures influence the mechanical response of discrete assemblies. One of our hypotheses is that a more physically based calibration of a material may ultimately result in a more unique set of contact parameters. We performed numerical experiments using (1) varying clusters made up of a finite number of particles to represent unique shapes, and (2) using two different assembly generation algorithms to pack them. These issues of the competing roles of particle packing and contact parameters become even more important when modeling coupled processes such as fluid flow and mechanical behavior. This is because of the obvious dependence of fluid flow on the material packing geometry.

The discrete element method has emerged as a powerful numerical simulation tool for a wide variety of industrial problems characterised by discrete/discontinuous features, in which different shaped objects with different sizes have to be considered. The pre-step in the discrete element simulation of practical problems often requires the generation of discrete objects packed in a form that can represent various realistic situations, an issue which appears not well addressed elsewhere. In many practical situations, disks/spheres may be sufficient to represent actual situations but often different size disks/spheres have to be used in the simulation. In other situations, more diverse shaped objects such as ellipses and polygons are necessary. The random packing of disks/spheres is an issue that has attracted a reasonable amount of attention over the last few decades, but nevertheless, preparation of an initial distribution of a large number of such discrete objects in a realistic (random) manner is not trivial. The existing packing algorithms are mainly of a non-constructive nature and therefore may involve substantial computer costs for large scale situations. Furthermore, very limited work has been reported with regard to the development of efficient packing algorithms for other shaped objects such as ellipses and convex polygons. In Feng et al., a novel numerical method is proposed to constructively generate a realistic random packing for different sizes of disks within a domain, motivated by the idea of the advancing front technique employed in mesh generation procedures. There are two different implementation versions of the algorithm, termed the closed form and open form. The closed form ignores the existence of the boundaries of the domain or container to be filled, while the open form takes the restriction of the boundaries into consideration when generating the disks. The algorithm has proved to be very effective, having a linear complexity and taking only several seconds to generate 1,000,000 disks with graded size distribution. The main purpose of the present paper is to extend the same methodology to other shaped discrete objects including elliptic particles and convex polygons, and as well to 3D spheres with different sizes. Computational issues unique to each special case are addressed briefly.

Over the past three decades, discrete element modeling has become an increasingly popular method to analyze a wide array of discontinuous systems. Since its genesis, the discrete element method has been widely used in the geotechnical field for describing the behavior of particle interactions. Much of the research in discrete element modeling has since focused on the development of faster and more computationally efficient algorithms to accomplish DEM simulations with the goal of modeling physically significant systems. However, considerably less research has been devoted to validating the microscopic behavior of the model in relation to physical systems. Cundall and Strack originally proposed using circular disks and Newtonian interactions to simulate the behavior of particles. The use of disks to approximate granular behavior continues to be a popular analysis method in the field. However, perfects disks are stable in a hexagonal packing configuration, as shown in Figure 1. A corresponding stable surface in a hexagonally packed configuration forms a 60° slope with respect to the horizontal. Recent work by Zhou et al. has shown agreement between empirical and numerical results for piling of spherical glass beads. However, this research does not address the influence of dense packing configurations. Most particles of practical engineering concern are not perfect disks and do not often exhibit hexagonal packing. Eliminating this behavior in the simulation environment is important in accurately representing many physical systems. Model behavior inconsistent with physical behavior at the microscopic level can propagate throughout the system and lead to overall incongruities between the simulation and the physical system it attempts to represent.

In most discrete element modeling, the preliminary effort focuses to create an object to be simulated, which consists of many particles. Unlike the finite element mesh, the initial configuration of particles to form an object cannot be generally determined only from the geometrical information. Except in some special cases such as a particle packing with very loose packing density and a packing arranged in regular lattice, an initial stable particle packing is achieved by (1) generation of particles with specific shape and size, and by (2) arrangement of the particles. Usually, the first procedure is performed from the computational geometry information to generate a pre-initial packing of particles. For the latter procedure, except some cases, there is still no analytical method or golden rule to determine the position and direction of each individual particle in a packing, especially for dense packing of particles with wide range of size and shape. Hence, it generally requires using mechanical ways, and consequently DEM itself is generally used. Depending on the pre-initial packing condition, this initial packing simulation sometimes consumes longer time than the total time of the simulation of main interest. On the other hand, the macroscopic material properties, such as stiffness and yield strength, in discrete element modeling are characterized not only by microscopic parameters but also by the configuration of the particle arrangement after initial packing. It is due to the space distribution of inter-particle contacts or bonds which form a force passing network. In other words, we cannot fully predict the bulk properties of many particles system from the information of single particle or of two or three particles before packing. It is, therefore, important to establish the relationship between packing methods and the properties of resultant packing for various cases in discrete element modeling in relation to the calibration of the model parameters. Furthermore, we need to investigate the mechanisms of packing behind the system of particles with certain size and shape.

The discrete element method is recognized as a powerful tool for studying granular materials. The first generation of discrete element models idealized particles with discs in 2-D and by spheres in 3-D. Later, polygons were used to improve particle shape idealisation. However, polygon elements are demanding on computational time. Ting et al. reported an increase in execution time of at least one order of magnitude for polygons compared to circular or disc shaped particles. Lately, the use of clusters of particles has been pursued by researchers. This approach does not require much modification of contact detection schemes usually used with circular shapes. However, the composition of analyzed samples in terms of percentage of each association is taken arbitrarily, which usually is unrepresentative of actual particle size distribution of granular materials. This paper presents the use of digital images to provide a real packing configuration for samples and then uses the cluster concept to improve particle modeling for use in discrete element analyses.

Virtual Reality (VR) based methodology has important potential applications in engineering. Among the earliest engineering applications of VR are the advanced VR simulators with realistic force feedback in civilian and military aircraft design and in pilot training. Rapid advances in computer hardware, software and information technology now enable the use of VR environments with realistic force feedback for design and prototyping in a wide range of engineering fields. In the area of vehicle design, VR applications include testing for crashworthiness, ground-vehicle interaction and driver visibility. An important class of potential applications of VR environments with realistic force feedback is in machine-medium interaction problems. This class of problems includes the machine-medium interaction in heavy earthmoving equipment and off road vehicles. VR environments with realistic force feedback are potentially useful tool in the design cycle for new equipment. Early attempts in this field have had limited applications only in examining the issues of driver field of view. The primary difficulty in these early applications has been the lack of a realistic and real-time modeling of the medium-tool interaction. Realistic mechanic based models of granular materials have existed for a number of years. Discrete element method has been successfully used to model the flow of granular materials. The models are capable of correctly representing the mechanical behavior of masses of granular materials. However, these models cannot be used in virtual reality environments since they are computational very slow. We are proposing a new approach to this problem by using neural networks. The methodology is based on the premise that neural networks can learn to accurately describe the instantaneous motion of the particles in amass of granular material. The methodology contain some DEM-based concepts and will be trained using data from DEM simulations.

The Java programming language and the Java3D graphics application programming interface (API) provide certain benefits over traditional programming languages and computer graphics libraries in the development of engineering applications, such as the implementation of discrete element methods (DEM). Java, as a pure object-oriented programming (OOP) language, enables the implementation of more robust and extendable programs through its information hiding, encapsulation, inheritance, and polymorphism characteristics. In addition, it provides portability, architecture-neutrality, and multithreading capabilities. This paper addresses performance concerns about the use of Java that arise due to the intensive computational demands of the DEM. Furthermore, the utilization of database technology to efficiently manage the huge amounts of data that are typically generated by DEM programs is discussed. An implementation of combined DEM and finite element methods (FEM), is presented as it utilizes the Java and database technology.

Many discrete element programs are not user friendly. They are usually run in batch mode that does not allow user interaction. Problems during an analysis simply result in abortion of the execution. Input procedures are usually scripts. Command syntax has to be memorized and results are seen only after program execution has been terminated. Even some commercial programs that have Windows-style pre- and post-processors are not easy to understand, and require a rather steep learning curve in their use. Input files can be saved, avoiding re-typing, however all input calculations, such as element generation, contact model calculation, etc., are repeated every time the program is started. Output typically consists of animation of element motion, shown simultaneously during computation. Additional analyses can be performed later after the calculation is finished. For an average technically educated user, the lack of click, drag-and-drop mouse activities is typically taken as a disadvantage and, in our opinion, presents one of the main obstacles for wider spreading use of DEM, not only in industrial application but also within the academic and research community. In an attempt to overcome this problem, a graphical user interface for a two-dimensional distinct element program DEMEOS has been developed. This program is developed to model the excavation of the Alberta oil sands, a rather peculiar viscous-frictional material with time and temperature dependent behaviour. DEMEOS is coded in Sun Java JDK1.3.1. Its functional flowchart is shown in Figure 1. Java was chosen because of its increased portability, but it has its pros and cons. It is a language under development, fast changing, and requires frequent re-coding from developers in order to keep pace with the current trends and capabilities. Comparison between various programming solutions in DEMEOS (especially the development of customised Java data structures) and an exemplary C++ DEM program like PFC will provide an interesting discussion in determining the advantages and disadvantages of one over the other.

The research here presented has employed the newly evolved combined finite-discrete element method in the development of novel numerical solutions for the analysis of failure and collapse of reinforced concrete structures under hazardous loads such as those arising from blast loading. Recent bomb explosions have made engineers more aware of the need to develop efficient computational tools as a viable cost effective means of analysing the failure, fracture and collapse of structures and structural elements. The first step to achieving this was the study of the structural response, failure, and collapse of individual structural elements. Thus the research in this area is taken further by using numerical solutions to study the behaviour of reinforced concrete beams to the point of failure. The results are implemented into the combined finite-discrete element method through a novel computationally efficient two noded beam element. Numerical integration across the cross section of the beam element is applied to facilitate the application of non-linear constitutive laws for both steel and concrete for the case of multi-axial bending coupled with axial force. The accuracy of this new element is tested and validated under both static and dynamic loading situations using analytical solutions together with experiments undertaken at the University of Alberta and The Swiss Federal Institute of Technology, Zurich. The proposed element has the advantage of reducing the size of the problem by fifty percent through elimination of the rotational degrees of freedom using static condensation. The new element, enables the same finite element mesh to be used for the discretised contact solutions, thus further reducing the CPU time required. When implemented into the combined finite-discrete element method, the proposed numerical solution also takes into account contact-impact and inertia effects.

In discrete element methods the simulated bodies are typically assumed to be infinitely rigid in order to reduce the computational cost. However, there are multibody systems where it is useful to take into account the deformability of the simulated bodies in order to enable the evaluation of their stress and strain distributions. This paper focuses on the simulation of systems of multiple deformable bodies using a combination of discrete and finite element methods (FEM), with some simplifying assumptions that are necessary to make the solution of the problem feasible. In traditional mixed FE formulations the contact effects can be taken into account using Lagrange multipliers methods and keeping the contact surfaces and forces as unknowns together with the unknown displacements. This approach results in huge systems of highly nonlinear coupled equations due to geometric as well as boundary nonlinearities. Furthermore, the parts of the bodies that may come in contact, typically, have to be identified before performing the simulation. However, no prior knowledge of the upcoming contacts is available in the multibody systems under consideration. Considering the excessive computational requirements, due to the huge number of degrees-of-freedom (DOF) and the high nonlinearities of the coupled systems of equations, it is unrealistic to solve problems involving many interacting bodies using such classical contact FE approaches. Simulations of deformable bodies with reasonable computational cost are enabled by incorporating FEM in DE analyses using certain assumptions that uncouple the contact interactions from the equations of dynamic equilibrium. In particular, the DEM are employed to identify, at each simulation step, the bodies in contact and determine the contact forces. Then, either a FE or a DE formulation is used at the individual body level to describe the equations of motion, depending whether the body under consideration is deformable or rigid, respectively. In case of a deformable body, the strains are assumed to be sufficiently small to permit a small strains, although large displacements, analysis. The deformability of individual bodies is taken into account using a displacement-based (DB) Updated-Lagrangian (UL) Finite Element (FE) formulation. Finally, an explicit time integration method, specifically the central difference method (CDM), is used to perform the numerical direct integration of the equations of motion and determine the new displacements, as well as the deformations and stresses wherever needed, of each body. Having computed the motion of each discrete body at a new time step, the positions of all discrete bodies are updated and a new contact detection process determines the new contacts and evaluates the corresponding contact forces, which are then used in the following time step.

This paper demonstrates that a mesh based partition of unity method, known as the manifold method, provides a unified framework for the analysis of both discrete and continuum problems. At the core of the method is construction of a mathematical mesh upon which finite covers are generated. Approximation functions are then developed upon these compact covers. The construction scheme applies for simple discrete element as well as continuum element that are even marked with discontinuities. The procedure as well as sample applications are presented.

The development of the finite element method (EEM) is recognized as an important contribution to numerical analyses. FEM has been used extensively in the structural analysis of continuum and other numerical simulations. During the past, numerical discontinuous models have been developed for contact impact problems. The "joint element" was developed by Goodman et al. to model the discontinuous joint interfaces, and the "distinct element method" developed by Cundall has been used for numerical computations of jointed or blocky rocks in many engineering problems. Based on the formulations of variational principles, contact mechanics have been proposed for contact finite elements in modeling discontinuities, such as the penalty method, the Lagrange multiplier method, the perturbed Lagrangian method and the augmented Lagrangian method. DDA considers both static and dynamic computations. Since the static condition is ultimately in the equilibrium state after a certain period, time steps are used for both cases. Dynamic calculation becomes natural as time is taken into account, and for static calculation, the process is assumed to be slow as the initial velocity of each time step is set to zero. The inertia matrix is more important than the elastic stiffness matrix in large displacement and deformation computations because the deformation strains are the first derivatives of the displacements, and the integration part of the inertia matrix is two orders greater than that of the elastic stiffness matrix. An accurate integration of inertia matrix allows the description of complete wave-propagation equations and prevents the strain components from sudden accelerations in every direction. Since the strains are well controlled, the process of adding or removing springs can only disturb the block system locally around contact positions. As a result, the open-close iterations converge easily and the multi-body discontinuous computation is stabilized.

The fracture behavior of solids such as concrete and rock under impact loading has become an area of increased research in recent years. These materials are complex and extremely heterogeneous, especially after they degrade from solid to particulate. Continuum based methods such as the finite element method are challenging to apply to this class of problems, and are plagued by the need for continuum constitutive models, severe element distortion and frequent remeshing. Furthermore, erosion of elements implies loss of mass, which may alter the solution for problems with extensive damage. This paper presents a discrete element method (DEM) for the analysis of such problems, in which solids are allowed to progressively damage into particulates. The study is conducted using a modified version of the program TRUBAL, originally developed by Cundall and Strack, for modeling the mechanical behavior of granular materials. Jensen et al. introduced the concept of clusters to obtain a better representation of particle shape for modeling granular materials. Clusters are formed by combining a number of circular-shaped particles in a semi-rigid configuration. In addition to our work in clustering, there has also been significant new activity by others along these lines, such as Qiu and Cruse, Sawamoto et al., and Thomas and Bray.

Our understanding of particle-fluid dynamics has been severely limited by the nonexistence of a high-fidelity modeling capability. Continuum modeling approaches overlook the microscale particle-fluid interactions from which macroscopic system properties emerge, while experimental inquiries are plagued by high costs and limited resolution. One promising numerical alternative is to simulate particle-fluid systems at the grain-scale, fully resolving the interaction of individual solid particles with other solid particles and the surrounding fluid. Until recently, the direct simulation of these systems has proven computationally intractable. In this research, a robust modeling capability for the direct simulation of particle-fluid systems has been formulated and implemented. The coupled equations of motion governing both the fluid phase and the individual particles comprising the solid phase are solved using a highly efficient numerical scheme based on the discrete-element (DEM) and the lattice-Boltzmann (LB) methods. Particle forcing mechanisms represented in the model to at least the first order include static and dynamic fluid-induced forces, and intergranular forces from particle collisions, static contacts, and cementation. The coupled method has been implemented into a generalized modeling environment for the seamless definition, simulation, and analysis of two-dimensional particle-fluid physics. Extensive numerical testing of the model has demonstrated its accuracy over a wide range of dynamical regimes.

Discrete particle simulations (including DEM simulations) inherently have a critical limitation: the number of particles. Several attempts have been made to overcome this limitation. Kazari et al. proposed a similarity model for DEM simulations of fluidized beds, in which large model-particles are used instead of small physical-particles, and showed that the bubble motion in a 2-D fluidized bed is well simulated by their similarity model. Sakano et al. proposed an imaginary-particle model in which the assembly of small physical-particles is represented by large imaginary-particles. Takeda and Horio proposed a hybrid model of DEM and two-fluid model, where the inter-particle stress in the constitutive equation of two-fluid model is evaluated using a DEM simulation in the corresponding small sample grid. In these models, the similarity model by Kazari et al. seems to express well the flow behavior in gas-fluidized beds in spite of its simplicity. Kazari et al. pointed out that flow behavior is changed when the size of model-particle is too large. However, the physical reason of this deviation and the limitation of their model is not clarified. In the present study, the limitation of the similarity model by Kazari et al. is studied by applying it to a DEM simulation in a 2-D fluidized bed. The calculated bubble motion is validated with the experimental result. We investigate the effect of the size of calculation particles on the bubble motion.

Tokyo has a number of river embankments in below-sea level areas called "Tokyo lowlands." Because the deposits in Tokyo lowlands generally exhibit a thick soft soil layer, it is very important to study the liquefaction behavior for the mitigation of earthquake disaster. It is well known that the major cause of Yodo-gawa river dike damage by the 1995 Hyogoken-Numbu earthquake in Osaka, Japan, was the liquefaction of the ground where the loose sand layer existed. For these reasons, an earthquake-proof project of embankment is carried out by Tokyo Metropolitan Government presently. Figure 1 shows the severely damaged cross section located near the Yodo-gawa river mouth within about 40km from the epicenter. The concrete parapet wall at the top of embankment slipped into the riverbed and the maximum vertical permanent displacement was about 3m. Recently, the rapid increase of computer abilities has drastically extended availability of Discrete Element Method (DEM). In the field of geotechnical earthquake engineering, the liquefaction is one of the most important, complex, and necessary topics. In this study, the liquefaction analysis method by using DEM (Code name: PFC^2D) in which the excess pore water pressure is taken into account to estimate the large permanent displacement due to liquefaction at the microscope level, is applied to the Yodo-gawa river dike damaged by the 1995 Hyogoken-Numbu earthquake.

The application of DEM in fluidised bed modelling was first pioneered by Tsuji et al. and since then, most of the work in the literature utilising this tool have either been 2D or pseudo-3D (3D discrete particles, 2D continuum fluid) simulations. In this paper, fully 3D and 2D DEM-based discrete particle-continuum fluid models are used to simulate some aspects of uniform fluidisation of a Geldart group D particle bed and the results are compared.

A coupled method between discontinuum and continuum approaches has been developed to simulate granular particles moving in a flowing fluid. The movement of granular particles is modeled using the discrete element method, while the finite volume scheme is used to simulate the fluid flow. The existence of granules in the flow region is assumed to cause a reduction in the amount of fluid flow in that region. The flowing fluid is expected to induce drag forces on the granules which are computed based on the relative velocities between the fluid and the particles. Buoyant forces on a particle are equal to the weight of the fluid displaced by the particle. All forces acting on each individual particle are summed in determining the particle movement. An example is given to illustrate the possibility of the proposed method in simulating granular particles moving in a flowing fluid. The problem of modeling wet granular flow requires the solution of two major physical problems; granular material as discontinua and fluid as continua. The discrete element method or DEM, proposed by Cundall and Strack, conceives granular materials as the assemblage of distinct rigid particles. Particles can interact with each other or with a solid boundary only at contact points. By using the contact force-displacement law, the forces at the contact are related to the magnitude of the overlap between particles or between particle and a solid boundary. The displacements and velocities of each individual particle are calculated from the summation of all forces acting on the particle using Newton's Second Law of Motion. The movement histories of individual particles at each time-step are traced based on the force-displacement law and the Newton's Law. Various particle shapes can be used such as polygonal, elliptical, circular shape, etc. The circular shape is selected in this paper due to its simplicity. For the fluid computational part, fluid flow is modeled using a finite volume scheme. The flow field is divided into a finite number of subregions, called cell. The continuity and Navier-Stokes equations are both applied at each cell in the flow field. To couple the movement of particulate material into fluid flow, fluid/particle interactions must be taken into account. Cundall classified five types of fluid/particle interactions depending upon the situations. Babic' and Shen as well as Sun and Vinogradov considered drag forces as the interaction of fluid flow to solid particles. Chan simulated the wet granular flow by considering buoyant and drag forces as the interaction of fluid flow to particles; on the other hand, the interaction of particle to the fluid flow is represented by using the different permeability to control the flow volume. In this paper, drag and buoyant forces are included in the fluid/particle interaction. For interaction between granular particle and the flowing fluid, the existence of the solid particles will cause resistance to the flowing fluid, In essence, the discrete element method is used to track the particles movement and the location of the particles is taken into account in the continuum calculation of fluid flow by introducing local impedance. The continuum method obtains the position of granular particles, computes the flow velocity and pressure, and imposes a drag force on each particle for the DEM calculation.

It has long been recognized that powders show four distinct types of behaviour, i.e., cohesive, aeratable, sand-like and spoutable, when fluidised. These effects are believed to be due to the variation in the relative importance of interparticle cohesive forces, such as van der Waals, capillary and electrostatic, compared with the drag force exerted on particles by a fluidising gas. However, such effects are difficult to quantify experimentally due to difficulties in the force measurements. In recent years, the Combined Continuum and Discrete Model (CCDM) has been developed to study flow of particulates where not only interactions between particle and particle but also particle and fluid are important. In CCDM, the continuum gas flow is solved by local averaged Navier-Stokes equations and the discrete solid flow by Newton's second law of motion based on Discrete Element Method with their coupling by Newton's third law of motion. This paper presents a study of the effect of interparticle cohesive forces on gas fluidisation by means of CCDM.

Permeability is an important material property of porous and granular materials such as rocks, sands, and soils. The accurate measurement or estimations are very important in modeling the flow through porous media as its relevance to hydrology, soil mechanics, oil industry, and environment protections. There are two fundamental ways to obtain the estimation of permeability: experimental measurement and estimation via theoretical modeling. The experimental method measures the flow through a small sample in a laboratory, then the permeability is calculated. The permeability obtained this way is the averaged result from the representative volume and is not generally applicable in real large flow simulations. The upscaling is to try to bridge the gap between laboratory results and the field parameters. The theoretical estimation of permeability provides a good approximation by using an idealized networks of tubes/pipes through which the fluid flows, and the flow in the pipes is considered as 1-D flow with prescribed pressure or boundary conditions. The networks are generally treated as regular although recent studies begin focus on the fractal structure or statistical features. However, it is very difficult to model the detail flow patterns in random media where different sizes of particles and different porosities are presented. The common approach is to use an averaged representative volume with regular averaged particle size and uniform porosity. Clearly, this is far from the reality in the random porous media. This paper presents a new discrete element approach to the permeability modeling using spherical particles saturated in viscous fluids and treating the interactions among the particles and between fluid and particles in a discrete element manner. For a given particle size distribution, the porosity is calculated. By using the proper boundary conditions, the liquid flux is computed, and the permeability is then estimated. This will in turn give a relationship between permeability and porosity. The simulated results via discrete element method will be compared with experimental results or well-known Kozeny-Carman equations.

The dynamic flow of water and other viscous fluids through granular soils is commonly modeled using Darcy's law with strain-independent permeability coefficients. Experimental studies show that this law is not valid for large nonlaminar fluid velocities which may develop under high hydraulic gradients. Furthermore, particle rearrangements leading to substantial variations in porosity are prevalent during soil liquefaction and other extreme flow driven phenomena, such as piping. These variations in porosity affect significantly the soil hydraulic conductivity and stiffness characteristics. Experimental investigations furnish only macroscopic phenomenological characterization of these complex response mechanism. Microscopic numerical analyses provide effective alternative tools to assess the hydraulic and deformation characteristics of granular soils in a coupled fashion. Herein, pore water flow is idealized using averaged Navier-Stokes equations which are coupled with a discrete element model of the solid phase. Such coupled formulation is computationally tractable and provides adequate level of information on pore water flow. This paper presents preliminary results of an ongoing research effort to simulate the micro-mechanical fluid-particle response under extreme deformation and flow conditions such as liquefaction and piping.

Semi-autogenous and autogenous milling has become an integral component of modern mining operations. They allow large volumes of raw rock feed to be processed efficiently and cost effectively. Current semi-empirical methods for the design of mills are based largely on data obtained from pilot and full scale mills. While these methods are highly successful and currently indispensable, they provide little insight into the mechanics of the charge motion and limited understanding of the fundamental processes that occur within the mill. The Discrete Element Method (DEM) is a promising numerical tool for the simulation of the complex dynamic particle motion and interactions within mills. It is envisioned that DEM will eventually be used in conjunction with empirical methods to better optimise the mill design process. Prior to this occurring, however, DEM must be rigorously validated from its most fundamental level upwards. This paper presents such a validation, using three dimensional trajectory data from a pilot mill against which to benchmark DEM.

The behavior of granular packing is of fundamental importance to soil mechanics. In particular, there has been substantial interest in the stress distribution within a granular assemblage under shearing in dealing with its localization behavior. The numerical description for such a stress can be derived by the contact force distribution based on a discrete computation like DEM, since it is difficult to solve the equilibrium in the regime of large deformation with the continuum approach. For the experiments on the granules, e.g., trap-door tests are to be considered herein, a photo-elastic analysis is the only way to capture the profile of the contact forces. An inverse analysis strategy is an alternative way to understand the mechanical response of the packing. For example, the Equivalent Inclusion Method is used to estimate the stress field in the granular assemblage using the observed displacements of granules and the traction along the boundary under the trap-door tests using aluminum rods with different radii. This type of inverse analysis is advantageous for estimating the inelastic stress distribution without full knowledge of the constitutive equation of the material. The aim of the current paper is to examine the stress distribution of an assemblage of circular disks by both an inversion, using the observed displacement profile of the loci of the rods during the tests, and a conversion through the DEM computed contact forces among the rods. Measurements of the loci are taken at the center of the rods through the image processing system and on the traction by the deformable tip with a strain gauge. The DEM computation is made under the same conditions as the experiments.

Validation using high-quality physical tests is important to establish the ability of discrete element methods to capture accurately the response of granular materials. A series of triaxial and plane strain compression tests on steel spheres with face centered cubic and rhombic packings was performed. Simulations of these tests were then carried out using the distinct element method (DEM). While the numerical simulations captured the observed laboratory response well, the post-peak response is shown to be sensitive to the coefficient of friction assumed along the boundaries. Simulations of physical tests on hexagonally packed rods have been useful for providing insight into the response of two-dimensional granular materials. In this earlier study, biaxial compression tests were performed on uniformly sized rods and on rods with a known size distribution. Recognizing that real granular materials are three-dimensional, this work has been extended to three dimensions. The study considered here examined the response of uniform spheres with both a face-centered-cubic (FCC) packing and with a rhombic packing in compression tests. As with the two-dimensional case, the sensitivity of the response to minor variations in ball size distribution was examined. This paper firstly outlines the physical test procedure, the numerical simulations are then described, and finally some representative results are presented.

The arching phenomenon is manifested through the reduction of stresses experienced by underground structures. Arching plays an important role in geotechnical engineering construction, such as excavations, retaining structures, pile groups, tunnel boring machines, culverts and underground facilities. The arching mechanism is intrinsic to granular material/rock mass independent of scale effect. Its fundamental mechanism relates to the ability of discrete units to transfer loads through interaction in a preferable geometry, and thus to bridge between the zone (or point) of load application to the zone (or points) of reaction. Though investigated extensively, very little work has been done to address the fundamental mechanism related to the discrete particles interaction. Testing results of an advanced experimental technique is presented and analyzed. A model of granular material made of photoelastic particles is utilized. The model and the sophisticated image and global data acquisition system allow one to track the development of the arching within ideal granular material during a trap door experiment by following the motion of each particle and the contact forces between the particles. Visual and quantitative analyses are presented demonstrating the relationship between the global arching phenomenon and the particle interaction on the micro-level. The information allows one to observe the changes associated with the arching mechanism and the stress variation resulting from it. Comparisons are possible between DEM and the test results.

In recent years considerable effort has gone into applying distinct element methods to tumbling mills, and for good reason, as insights gained through DEM modeling promise significant improvements in both energy efficiency and general performance. These DEM codes use soft contact between particles and generally assume the operation of normal and shear springs for motion perpendicular and tangential to the contact plane, with the latter shear spring giving way to simple frictional sliding under excessive shear forces. The springs can be linear or nonlinear. In addition several damping models have been invoked, perhaps the simplest, borrowed from fluid mechanics, being velocity dependent damping introduced independently to each spring system, creating normal and shear spring-and-dashpot contact systems. However, several features commonly present in DEM codes and of appreciable value in DEM modeling of coherent solids would seem to be inappropriate for the modeling of tumbling mills.

Discrete element method is an ideal numerical technique for studying the constitutive behavior of not only granular materials, but also cohesive materials such as a clayey soil. Presented in this paper is a three-dimensional discrete element method for modeling the behavior of saturated clays. In addition to the mechanical interactions as between two granular particles, there are physico-chemical interactions between two clay particles. Two primary physico-chemical forces are the van der Waals attractive force and the double-layer repulsive force. Furthermore, unlike granular soils, which are bulky, clay particles are plate like. In the present study, a three-dimensional discrete element method is developed and used to conduct one-dimensional compression of an assembly of particles and some results are presented and discussed. Micromechanical methods such as the discrete element method (DEM) are very useful in studying the fundamentals of the stress-strain behavior of clayey soils. Proposed originally by Cundall and Strack to analyze a collection of granular particles, the method has since been used to study many other types of discrete, and in some cases continuum, materials. Anandarajah developed a two-dimensional discrete element method for the analysis of clays, where particles were approximated by rectangles. Procedures have been developed in recent years to quantify the physico-chemical forces between clay particles, which were used in the two dimensional studies. A more realistic simulation requires a three-dimensional analysis methodology, and this is achieved in the present paper.

Recent experiments by Eccleston and Miller, have directly measured strongly attractive particle-particle interactions by gluing 15 |j.m diameter Borosilicate glass spheres to the ends of cylindrical glass cantilevers and measuring the cantilevers' deflections as the particles are brought together and pulled apart. These experiments show a wide range of qualitative behavior of clusters of particles; the most curious of which is stick-slip, which would not be expected for force interactions that are between particles' centers of mass. We assume a van der Waals attractive force between macroscopic spheres and simulate the experiments using a Discrete Element Method DEM. It is found that the van der Waals force approximation accurately returns the values of the cantilever deflection measured in the experiment. We also find that the stick-slip behavior is recovered from the DEM model when the particles have a non-zero coefficient of friction.

Clay soils are commonly formed by the slow deposition of clay particles in water, e.g., lacustrine or marine bodies. Within these environments, the particles settle and eventually consolidate to form soil strata. During the settlement process the water chemistry and clay mineralogy will greatly influence the interaction between particles, i.e., whether particles will settle individually or in clusters known as "flocs". Initially the structure of the soil will reflect the settling behaviour, e.g., a flocculated soil will have a more open structure, with a greater volume of voids between particles, compared to a dispersed soil. As more particles accumulate, the stresses associated with the self-weight of the particles will increase and will begin to dominate the mechanical behaviour and structure of the material. Many of the clay platelets will eventually develop preferred orientations (with the longest axis orthogonal to the major principal stress direction), creating an anisotropic structure with variations of the mechanical behaviour in different planes. Despite this simple conceptual model, the contact points and stresses through the soil body are very complex and vary considerably during loading and unloading, and it is thought that particle clusters still dominate the behaviour and act as crushable, compressible "macro-particles". Regardless of this complexity, the soil behaviour is still controlled by the applied stress field and the materials alone, and this structure is generally referred to as a sedimented state. Much of our current understanding of the mechanical behaviour of clay soils comes from laboratory tests on reconstituted soils with this type of structure; hence this is also termed the intrinsic state.

Self-assembly is a promising process from fundamental and industrial points of view. A bottom-up type process for the fabrication of nanostructure is much less expensive and much more efficient in energy consumption than a top-down type process, such as photolithography. For example, two-dimensional arrays of colloidal particles can be fabricated by drying process of colloidal dispersion. These materials are applicable to optical devices, lithographic masks, and so on. However, it is still difficult to obtain highly ordered arrays in large area by drying process. To improve the degree of the order first, one needs to understand the dynamics of evaporation-induced self-assembling process dominated by capillary immersion force. We modeled self-assembling process of colloidal particles induced by capillary immersion force with Discrete Element Method. This method has been applied to various phenomena related to powder technology. The concentration of colloidal dispersion and the frictional force between colloidal particles are varied in this study. The dependence of the degree of the order for two-dimensional array is studied using Voronoi Polygon analysis.

The failure mechanisms of factional granular materials are characterized by complex failure modes and, furthermore, they show a highly anisotropic bias due to their inhomogeneous microstructure. Classical continuum mechanics formulations applying, for example, standard finite elements can consider only a small number of discontinuities and cannot encompass the entire fracturing process. However more enhanced continuum models, for example based on the Cosserat theory, may provide the possibility to account for a granular microstructure in a straight forward way. On the other hand, the mechanical behavior of granular materials can be well described on a smaller scale via discrete element models (DEM) if the motion of the individual grains is taken into account. The time evolution of such a model is followed solving the Newtonian equations of motion of the individual elements numerically. Our concern in this field of research is the relation between the microscopic quantities, like the contact forces and displacements, and corresponding macroscopic quantities, like stresses and strains via corresponding homogenization techniques. In this contribution we will propose a new technique and validate it by two-dimensional DEM simulations with polygonal particles. The evaluation of different boundary value problems proves the efficiency of the homogenization technique.

In this paper we investigate granular materials with respect to their response to shear. Certain number of particles are connected into straight line chains to create spatially extended, elongated grains. The effective friction coefficient of the material is measured as a function of shear speed and normal vertical stress. Non-spherical, elongated grains are investigated using a recently developed numerical technique for efficiently evaluating the Lagrange multipliers necessary for dynamics of particles with geometrically linear constraints. Material response to shear of granular systems with individual grains consisting of up to five particles is investigated.

In 1990, Babic et al. published a paper describing the stress and strain-rate relations in a simple shear flow of granular materials. They found that as the solid concentration increases keeping all other parameters constant, the stress at first was nearly proportional to the square of the strain-rate, then went through a sharp transition and finally became independent of the strain-rate at high solid concentrations. In 1991, Hopkins and Louge published a paper describing the microstructure of granular shear flows. It was observed that microstructures formed in shear flows. These microstructures were more prominent for inelastic particles. In 2000, Zhang and Rauenzahn published a paper describing the contact duration between colliding particles. They showed the effect of contact time in stress relaxation and demonstrated viscoplastic behavior of granular materials. Based on these findings it is conceptualized that stress and strain-rate relations depend on the internal structures of the granular materials. That is, stress and strain-rate relations for granular materials must include the evolution of the internal structure as an integral part of the constitutive relation.

Many interesting properties in a dense granular system, such as constitutive relations and effective heat conductivities, are not only determined by the particle volume fraction, but also the spatial and temporal correlations of the system. Spatial correlations in a granular medium as manifested by shear bands, force chains and size segregations have been observed and studied by many researchers. Effects of temporal correlations in a granular system have not received as much attention. Recent studies by Zhang and Rauenzahn suggest that the history of particle interactions is critical in understanding many interesting phenomena unique to slow and dense granular flows. These phenomena include viscoplastic characteristics, history-dependent stress strain relations, and phase transitions. To study the effects of particle interaction history, we built a discrete element simulation code capable of tracing histories of interacting particle pairs. In the numerical simulation, particles are placed into a cubic cell and then are subjected to random motions to create an initial configuration. Initial velocities of particles are specified according to prescribed shear motion. In the numerical simulations described in this paper, all particles are identical and spherical, the force model is similar to the one using by Jenkins and Strack. The normal force between particles is modeled as a parallel connection of a linear spring with coefficient K_n and a dashpot with resistance R. The tangential force is proportional to the tangential deformation. The magnitude of the tangential force is limited by Coulomb friction force determined by the product of the friction coefficient μ and the normal force. In the numerical simulation, the time is non-dimensionalized by 200m / K_n, where m is the mass of the particles. The stress is non-dimensionalized by K_n/(200a), where a is the radius of particles. The restitution coefficient for binary collision for the particles in the simulation is 0.8. The numerical scheme used in the discrete element simulation is described in recent papers by Zhang and Rauenzahn, and is not detailed here. To minimize the effects of boundaries, periodic boundary conditions are used in the simulation. In this paper, we briefly introduce the theoretical development of Zhang and Rauenzahn. One of the findings in the theoretical development is that after the granular medium is sheared for a sufficiently long duration, the age of a contacting particle pair obeys an exponential probability distribution for long contact ages. The relaxation of the contact stress in a dense granular medium is closely related to the relaxation time of this probability distribution. A constitutive relation derived from the study of the temporal correlation has shown many important features of granular material. In this paper, the author carried out several numerical experiments to examine these predictions by the constitutive relation. It is shown that the rheological behavior of a dense granular medium is vastly different from that of a polymeric material although both of them have finite relaxation times.

Crushable agglomerates can be simulated in the distinct element method (DEM) by bonding elementary spheres in `crystallographic' arrays. The behaviour of an element comprising these numerically-generated soil grains resembles that of a crushable anisotropic sand, in terms of irrecoverable compression, yielding and plastic hardening behaviour. Stress path tests are simulated to capture the response of the simulated soil to various combinations of deviatoric and spherical stress changes, δq and δp, in both loading and unloading. The plastic deformation is comparable with certain features of the Cambridge soil plasticity models, although the definition of a single yield surface separating elastic from plastic behaviour is arbitrary, exactly as it is with actual triaxial tests on real soils. Normality of the plastic strain increment vector to conventionally defined yield surfaces is not observed. It is shown that grain breakage is a crucial component of plastic soil behaviour. Yield surfaces obtained from numerical stress-path tests are best quantified by the percentage of bond-breaking (pbb) contours. These surfaces are compared with contours of yield defined alternatively from inspecting stress-strain curves. An "elastic wall" in voids ratio log mean effective stress (e, ln p') space is found, but the normality rule is not followed at yield. Similarity is observed in the shapes of yield surfaces at different degrees of strain-hardening, both in the e-ln p' space and the q-p' space. Their shape resembles the Modified Cam Clay model more than the Original. Lastly, a reduction of internal coefficient of friction, represented by the Cam Clay parameter M, is observed when stress-level is increased; this also resembles real soils.

Although the mining and quarrying industries have been using explosives for decades in breaking rock to allow the extraction of minerals, there is still much interest in finding ways to control the use of explosive energy to ensure each application is safe and the utilization of the loaded explosive energy is maximized. Stemming is known to be an effective control as a confinement to prevent the loss of explosive energy into the air through the collar before completing the work done in breaking and heaving the blasted material. An effective stemming should be able to retain an explosion gas pressure on the order of GPa. Stemming is a granular material made from sand or crushed rocks of different shapes and sizes usually passing a one-inch mesh. This paper demonstrates the use of a two-dimensional discrete element method to model the basic characteristics such as packing or density of a stemming material in terms of particle shapes and size distribution. A preliminary setup for the study of the interlocking capability in counteracting the stemming ejection pressure during blasting is also demonstrated. Described in Figure 1 are some sample cluster-shaped bodies of various sizes designed to simulate different shapes of stemming particles. More particles generated from a quadratic equation will be covered in the main paper. The acting forces within the stemming particles as they are set in motion are described below, in Figure 2, in terms of normal elastic and damping forces with adjustable spring and damping coefficients k_n and c_n respectively.

In the present study, the behavior of random assemblies of two-dimensional cohesionless granules of oval-shape particles under stress-control biaxial compression tests is investigated. In conventional discrete element method (DEM) proposed by Cundal and Strack, granular materials is considered as a system of a finite number of particles, in which equations of motion for both translation and rotation of each particle are derived based on the second law of Newton. Applicability of DEM to large deformation and localization type analysis has made it quite attractive and to date it is considered as one of the most reliable methods for studying the mechanical behavior of granular materials. However, huge amount of computational effort has been always considered as the major drawback associated with this method. As pointed out by Bardet and Proubet, the method uses explicit time integration, which is conditionally stable and converge very slowly. Moreover, an appropriate estimate for ideal damping parameters is required. The computational time can be quite cumbersome especially if particles with complex geometry are involved. In static or quasi-static problems where system parameters (such as location of the grains, and the magnitude and the direction of contact forces) need not be represented as a function of time, usage of dynamic equations and damping parameters, as proposed in conventional DEM, seems to be unnecessary. SDEM, originally proposed by Kishino for analysis of assemblies of two-dimensional circular disks, is an attempt to eliminate these unnecessary factors. The major advantage of SDEM lies in its high rate of convergence. Its major drawback is incapability of modeling dynamic problems. The methodology was later extended by Shodja and Nezmai for two-dimensional elliptical particles and employed in a series of strain-control biaxial compression simulations.

The compaction of powders in a cylindrical die is a typical method to produce molded products of powders and widely used in chemical, pharmaceutical, ceramic industries and so on. This type of compression is considered as a basic system of compression and therefore has been studied frequently by means of various experimental and theoretical methods by a number of researchers. Nevertheless, even now the knowledge and understanding of compression mechanism are not enough and difficult problems exist. For example, it is difficult to predict the stress and density distributions in the particles under the compression process. Further study is therefore needed. On the other hand, discrete element method (DEM) was proposed by Cundall and Strack, and has been widely used to analyze the discrete behavior of particles. DEM makes it possible to study particle-level effects on the behavior of particles such as velocity and stress distributions during dynamic processes. DEM have been used to analyze compaction processes of particles. However, there are some additional physical factors that can affect the behavior of particles in the compaction process such as roughness of the particle's surface. In this work, three-dimensional DEM simulations of compression test of particles in a cylindrical die are performed. The calculated results are compared with the corresponding experimental results. The comparison shows that the present DEM model, which implies smooth-surfaced particles, overestimates the stiffness of the particles' column. As the reason of this quantitative disagreement, the effect of surface roughness is considered. We propose a rough-surface model based on an observation by using an atomic force microscope (AFM) and applied it to the DEM simulation of compression test.

Powder compaction is widely used in the pharmaceutical, ceramic and automotive industry to produce near net shape components with complex geometry and high strength. Compaction induces very complex states of stress in the powder compact, which sometimes cause processing problems in industry and lead to undesirable non-uniform density distributions and even cracks within the pressed parts. Many experimental and theoretical studies based on continuum mechanics have been carried out to study the macroscopic mechanical responses of compacts. In recent years, the finite element method (FEM) has also been adopted in simulations of powder compaction. Although continuum methods (including FEM) are effective in studying the macroscopic densification process in compaction, it should be noted that the accuracy of the results strongly depends upon the faithfulness of the constitutive model and the quality of the material parameters. Furthermore, in continuum models, the microscopic properties of particles and the interactions between particles, which strongly influence the macroscopic behaviour, cannot be taken account. To get optimal mechanical properties of green compacts, better knowledge of the relation between powder characteristics and mechanical behaviour of the material during compaction is required. The Discrete Element Method, developed by Cundall and Strack, is an effective numerical tool to investigate the micro mechanics of granular materials. In DEM, each individual particle of an assembly is modelled separately and its motion is defined from the interactions with neighbouring particles. The detailed movements of particles can be traced. As a consequence, microscopic evolution of the internal stress and the structure of the assembly can be obtained. In the meantime, the macro behaviour of the whole assembly can be obtained by statistical means. In this paper, A 3D DEM code TRUBAL, which is capable of modelling friction, cohesion and local plastic deformation at the inter-particle contacts, is used to simulate quasi-static uniaxial compaction of a glass ballotini assembly. First, the critical time step that is proper for quasi-static simulations is discussed. The evolution of microscopic properties of the particle assembly, such as coordination number and plastic contact number, are investigated along with correlated statistical macroscopic behaviours such as stress-strain relations and anisotropy of the compacts. The link between particle interactions and mechanical properties of the compacts is thereby identified through the DEM results.

An experimental and numerical study was conducted on the geomechanical behavior of oil sands during penetration tests in which a steel wedge was forced into compacted oil sands. The boundary conditions simulated in the laboratory test are similar to many common excavation activities that occur during oil sand mining (e.g., shovel penetration into oil sands). The laboratory data was needed to improve calibration of input parameters for a two-dimensional numerical model of oil sands constructed with Particle Flow Code (PFC). The objective of the laboratory tests was to measure the force required to push a steel wedge into compacted oil sand. A wedge penetration test differs from conventional triaxial tests in that much more shear is generated. The measured force-displacement behavior of the steel wedge was then compared with PFC modeling results.

Many experimental studies have been undertaken to investigate the mechanics of the interface between granular materials and structural surfaces. The direct shear and shaft pullout tests are the most common. A variety of scales, configurations and boundary and loading conditions were used in these tests to capture the elements of the behavior. However, the majority of these investigations rely on measurements taken outside of the structure-medium interface to infer behavior at the interface. Effort has been made to correlate the macro-mechanical response of structure-medium systems to the microstructural characteristics of the interface. Uesugi and Kishida studied interface grain displacement using a reflex camera through a glass window of a simple shear device. Later, Boulon, working with a modified plane direct shear test and using stereophotogrammetic analysis technique, was able to track the movement of particles. Obtaining microstructural information using laboratory testing techniques is an inherently difficult, tedious and time-consuming task. The use of the discrete element method (DEM) to study and delineate phenomena dictating the behavior of the interface provides a valuable alternative since it takes into account the discrete nature of materials and gives detailed information at the particle level where interface behavior occurs. This paper presents some parametric studies on the effect of modeling parameters on formation of shear zones in structure-medium interfaces.

Many researchers have modeled the failure behavior of rock around underground openings. One approach is to assume an elastic-plastic continuum behavior for rock. To simulate stress-induced degradation, the material description should include softening plasticity. On the other hand, the numerical modeling of a softening material suffers from mesh sensitivity. To resolve this difficulty, the idea of non-local continua has been used by some researchers. In non-local plasticity, the stress at any point is a function of the average of strain around that point. Zietlow and Labuz have shown that the size of the averaging zone depends on the grain size. The distinct element approach is another alternative to model the softening i: behavior of rock. In the PFC^2D program, which will be used in this paper, the rock is considered as a conglomerate of rigid circular particles. The interaction between particles is elastic through normal and shear springs. If the normal or shear traction between two particles exceeds the bond strength, the corresponding contact will break. A broken contact does not have tensile strength, but in shear it behaves as a Mohr-Coulomb frictional material. By simulating the rock behavior as an assembly of distinct particles, crack propagation and softening behavior of rock around a circular opening in Berea sandstone will be studied. The results of the numerical test will be compared with those of a laboratory biaxial experiment.

In a previous study, two different discrete element methods (UDEC and DDA) were used to characterize the behavior of a blocky rock mass containing a rectangular-shaped underground opening. In eleven of the fifteen cases studied, both methods produced identical or very similar results. This paper presents the results of a follow-up study focussed on investigating the four cases in which the methods predicted drastically different behavior of the rock mass.

We are using DEM simulations in conjunction with experimental measurements of artificially cemented laboratory specimens to investigate the effects of small amounts of inter-granular cement on mechanical behavior. With current PC computing capabilities, it is possible to create DEM sample assemblies that faithfully replicate full three-dimensional laboratory test specimens. The ultimate objective of these studies is to devise procedures to predict the degree of cementation and strength of granular material with applications to wellbore stability, sand production from wellbores and liquefaction. Preliminary work has focused on wave propagation in sands and weakly cemented sandstones using Itasca's PFC Discrete Element Method (DEM) code. We discuss results of our DEM simulations concerning: (1) relationships between micro-mechanical processes and bulk deformation behavior, and (2) stress-induced velocity anisotropy around wellbores.

A numerical model for rock is proposed in which the rock is represented by a dense packing of non-uniform sized circular or spherical particles that are bonded together at their contact points and whose mechanical behavior is simulated by the distinct-element method using the two- and three-dimensional discontinuum programs PFC2D and PFC3D. The fundamental microproperties consist of stiffness and strength parameters for the particles and the bonds. Damage is represented explicitly as broken bonds, which form and coalesce into macroscopic fractures when load is applied. The model reproduces many features of rock behavior including elasticity, fracturing, acoustic emission, damage accumulation producing material anisotropy, hysteresis, dilation, post-peak softening and strength increase with confinement. These behaviors are emergent properties of the model that arise from a relatively simple set of microproperties. The fundamental particle in the PFC codes is a circle and a sphere, but complex grain-like shapes can be produced by bonding particles together and modifying the inter- and intra-grain properties. The macroscopic behavior of such models is affected not only by the deformability and strength of the grains and cement but also by the grain size, shape and packing. In this paper, we limit our study to arbitrary packings of circular grains joined at their contact points by a brittle elastic cement, and investigate the relationship between the grain/cement properties and the macroscopic properties measured from biaxial and Brazilian, tests on synthetic material assigned microproperties to correspond with Lac du Bonnet granite. Studies in which elongated bands of weaker and more compliant material were used to represent a schistous microstructure and thereby produce an anisotropic modulus and strength in both two- and three-dimensional models are described in Potyondy and Autio and Autio et al.

Rapid setting, thin, spray-on polymeric liner materials are being developed for underground rock support in Canadian mines. These materials have performance characteristics that lie between those of shotcrete and mesh. Two-dimensional numerical models using Particle Flow Code (PFC) are used to model thin spray-on liners in simulated laboratory tests. The objectives were to determine how to model liners with PFC and based on these results to gain some insight about the support mechanisms provided by a liner. Direct tension tests and block punching tests were numerically simulated and used to calibrate the PFC micro-mechanical input parameters. Based on previous laboratory testing, the representative properties for a generic liner material are a tensile strength of 2MPa achieved at 50% elongation. The tensile strength is equivalent to a rupture load of 2kN per millimeter of liner per linear meter. A typical thin spray-on liner thickness used to support rock around an excavation is about 4mm. Therefore, the PFC liner models were calibrated to break at around 8kN per linear meter of liner material.

Landslides are an important class of natural disaster that can lead to significant loss of life and significant property damage. Understanding when and where landslides may occur and under what circumstances is a key part of risk assessment and protective strategy formation. Understanding the path that any specific landslide may take and how far it will run from the initiation point are also important considerations. One class of large landslides is particularly dangerous because of the long distances that they can travel (more than several kilometres is quite possible) while exhibiting very low effective friction coefficients. Cleary and Campbell showed that this was possible (using a very simple two dimensional DEM model with periodic boundaries) because the shear in the landslide was restricted to a narrow layer at the base with most of the landslide mass travelling undisturbed over the top. This resulted in lower energy dissipation and the low effective friction that was required to explain this long run-out. Campbell et al. showed, using large scale DEM simulations of the mountain slope and valley flow (again in 2D), that the phenomena was more complex with the basal shear region being more blurred and continuous, but still restricted to the lower portions of the landslide. Here we extend the modelling to three dimensions and use some real topography to explore the effect of the additional realism on the DEM predictions of landslides.

Particle scale simulation of industrial particle flows using DEM (Discrete Element Method) offers the opportunity for better understanding the flow dynamics leading to improvements in equipment design and operation that can potentially lead to large increases in equipment and process efficiency, throughput or product quality. A range of large scale important industrial applications are presented here: 1) Assessment of sample bias for cross-stream samplers, 2) Segregation in a splitter bin, 3) Separation efficiency by a panel of a vibrating screen, 4) A grinding table for cement milling, and 5) Hard ore grinding: pilot and full size SAG mills. The DEM method used here is described in Cleary.

For the last eight years, Conoco Inc. has been funding a group at the University of Oslo to develop discrete element models (DEM) for the solution of a number of problems in the petroleum industry. A brief overview of the model and areas of application that we have found useful will be presented. Initially our interest was to predict faulting below the resolution of seismic data in order to better define the extent of faulted oil and gas reservoirs. We have found, however, that there are also applications for DEM models in improved seismic interpretation, estimating local stress variations for fracture density estimates and fracture stimulation treatments, and modeling well bore failures due to compaction and subsidence.

The ice pack covering northern seas is composed of a mixture of thick ridged and rafted ice, undeformed ice, and open water. Existing ice-ocean models of the Arctic ice pack are large-scale Eulerian continuum models that use a plastic yield surface to characterize the constitutive behavior of the pack. An alternative is to adopt a discontinuous Lagrangian approach and explicitly model ice parcels and the interactions between them. I have developed a granular model of the central Arctic ice pack that consists of tens of thousands of individual floes. Neighboring floes can freeze together to form viscous-elastic joints. Wind stress, temperature fields, and Coriolis accelerations drive the model pack. Simulations begin with all floes in the model pack frozen together and at rest. As the floes accelerate in response to the wind, fractures are initiated throughout the model pack and coalesce to form a fracture pattern that covers the basin. The fracture pattern defines many large plates 10–100 km in width that are composed of aggregates of the basic, small floes. Subsequent deformation takes place along the boundaries of the aggregates. Regional variations in the size of the aggregates correspond to regional variations in the gradient of the wind field and the thickness distribution. Once the joint between two floes has broken they can separate to form leads or be pushed together to form pressure ridges. Pressure ridging is incorporated via a parameterization based on discrete element simulations of the ridging process. At the end of each 24-hour period the model pack is reinitialized. Reinitialization begins by placing an undeformed configuration over the deformed pack. The thickness distribution in the deformed pack is passed to the undeformed pack New ice is grown in leads. The undeformed pack is run for one hour to create a new aggregate fracture pattern and to equilibrate the stress and velocity fields to the new daily wind field. The model is used to run simulations of several months duration. The simulation results are compared to the maps of Linear Kinematic Features (LKFs) of divergence, shear, and vorticity produced by the Radarsat Geophysical Processor System (RGPS). The output of the model pack is run through a virtual RGPS to generate similar maps of LKFs to use in the comparison. Regional variation of the basic floe size permits high resolution (~1 km) simulation of areas of interest.

Investigation of the granular flow induced by a single plough in a lab scale plough share mixer using PEPT has showed the development of two loops of circulation one on either side of the plough. The particle displacements are controlled by the number of blade passes. Comparison with DEM results showed good agreement with experimental results.

Particulate solids are ubiquitous in our daily lives. Uniform mixing and controlled segregation of particulate solids is the key to success of many industrial processes. However, the complex nature of particle interactions has defied a good understanding of particle mixing. A simple system such as a sand pile at rest with a slope less than the angle of repose behaves like solid material that remains at rest. If the same pile is tilted several degrees above the angle of repose particles start to flow like a fluid. This simple observation that embodies the richness and complexity of particle motion have led to years of research effort both from the theoretical as well as practical standpoint. While the physicists are trying to establish connections between granular dynamics and processes taking place on an astrophysical scale, many industries on a daily basis simply struggle with solid mixing to obtain products of an acceptable quality. The mixing process is typically carried out in relatively simpler equipment such as a rotating drum. Here mixing is influenced profoundly by the flow characteristics of the particulate matter to be mixed. The purpose of the paper is to study mixing and segregation of particles by elucidating some of the flow characteristics in rotating drum.

We present results from a parameter study investigating the stability of underground structures in response to ground shock. Direct simulation requires detailed knowledge of both the facility itself and the surrounding geology. In practice, however, key details of the geology, rock properties and reinforcement may not be available. Thus, in order to place bounds upon the predicted behavior of a given facility, an extensive series of simulations representing different realizations may be required. We have performed simulations of subterranean facilities using the Livermore Distinct Element Code (LDEC). This paper presents a description of the method, with emphasis on techniques for achieving improved computational efficiency, including the handling of contact detection and approaches to parallelization. In addition, some continuum approaches to the simulation of underground facilities are discussed along with results from underground explosions. Finally, results of LDEC simulations of dynamic loading of generic subterranean facilities are presented, demonstrating its suitability for this application. We employ the Distinct Element Method (DEM), as defined by Cundall and Hart. By nature, the Distinct Element Method can readily handle large deformation on the joints. In addition, the method detects all new contacts between blocks resulting from relative block motion. The Lagrangian nature of the DEM simplifies tracking of material properties as blocks of material move. It is also possible to guarantee exact conservation of linear and angular momentum. The joint models can be very flexible and can incorporate experimentally observed effects, such as cohesion, joint dilation, friction angle, and hysteresis. The DEM has been applied to a wide range of problems in geomechanics. For example, Antonellini and Pollard simulated the formation of shear bands in sandstone using the DEM. Morgan applied the DEM to the mechanics of granular shear zones. Heuzé et al. used the DEM to analyze explosions in hard rock. Cundall reviews the application of the DEM to simulation of granular material and rock.

In "wet granulation," liquid binder is sprayed onto the surface of a mixed granular bed, resulting in the formation of bonds between adjacent particles and thus the growth of aggregates. Such processes are common in a variety of industries, including the production of tablets in the pharmaceutical industry, the manufacture of fertilizer pellets, and agglomeration of wood fibres in the paper industry. The equipment used to facilitate wet granulation varies between industries and applications, but can be broadly classified as tumbling agglomerators (e.g., dish and pan pelletizers) or mixer agglomerators (e.g., paddle mixers, plough-share mixers and peg mixers). Common to all equipment is the need for improved means for controlling final aggregate properties (e.g., particle size distributions, porosity and aggregate strength). In general, these properties will be determined by the nature of inter-particle interactions and the distribution of liquid binder. Since these factors are essentially determined by the nature of the flows occurring during a given granulation process, it is reasonable to expect that improved control of granulation can be achieved through manipulation of parameters that effect granular motions. DEM can provide valuable information about such granulation processes, filling gaps in our understanding due to the limitations of experimental measurement. To begin, we perform dry (binderless) case studies of three well-known granulators to establish the overall flow patterns and to demonstrate the manner in which data provided by DEM can be used to identify which mechanisms control aggregation.

Bulk solids or granular materials, as they are sometimes called, can be defined as any material composed of many individual solid particles, irrespective of particle size and shape. Thus the term granular material embraces a wide variety of particle dimensions and shapes, ranging from the mixtures of ore and waste rock after blasting to the finest icing sugar. The first significant studies related to the storage and flow of bulk solids were reported at the end of nineteenth century. That work originated from the need to store large quantities of grain, and was concerned mainly with wall pressure affecting the structural design of silos and bins. The most famous researcher of that period was Janssen, who developed the formula for prediction of the wall pressure. His method is based on the concept of differential slices of infinitesimal thickness, and finite cross-section and perimeter. Lack of consideration of movement, shape, and blockiness of bulk materials has limited the application of Janssen's method to approximate calculation of wall loads on bins and bunkers, under static conditions. Even to date most design manuals present a set of empirically derived correction factors for use in conjunction with load predictions. In the past three decades, gravity flow of bulk materials has been studied by variety of methods, including: 1) Experiments, by analogy with the flow of other materials (e.g., sand) in bins and bunkers, 2) Full-scale field studies, 3) Large-scale physical models, and 4) Numerical modeling. The introduction of fast computers with large memory provides an opportunity for researchers to use numerical modeling to study the flow of bulk materials. Traditionally, continuum theory (Finite Element Method, FEM) has been used to describe granular material behavior. But many bulk-handling systems involve very large-scale, discontinuous dynamic or static behavior. Gradually, Discrete Element Method (DEM) with its unique ability for simulating the complete dynamic behavior of discrete, interacting bodies has introduced itself as a very powerful tool for numerical modeling of discontinua. DEM's applications in the area of flow of bulk solids have been extended from mining to pharmaceutical industries.

The mining industry routinely uses ball mills and semi-autogenous mills for processing ore bodies. A typical mine site processes 80,000 tons/day of ore employing two semi-autogenous mills and four-ball mills. Since the energy consumption is 2–5 kWh/ton of ore for grinding rock particles to sieve mesh sizes, a substantial portion of the cost of ore processing is expended in tumbling mills. The energy efficiency of tumbling mills can be examined by directly looking at the motion of rocks and steel balls inside the mill. The make-up of the charge and the lifter bars attached to the inside of the mill shell can be designed particularly to maximize the mass of ore fractured per unit of energy spent on the mill. At the same time, the unnecessary collisions of steel balls against the mill shell can be reduced. Furthermore, the cascading charge flow can be altered in such away as to maximize grinding efficiency. The emergence of the discrete element method (DEM) allowed simulation of charge motion in tumbling mill. In the last eight years, the DEM for simulation of tumbling mills has advanced sufficiently that it has become a very practical tool in the mining industry. This manuscript gives an overview of the DEM as applied to the tumbling mill problem. Both the two- and three-dimensional (2-d and 3-d, respectively) models as well as the parallelization of the 3-d code are described, followed by validation against plant size mill data.

Bulk solids handling is seen in many industrial fields. The dynamics of granular materials in industrial equipment or machinery, such as bins, hoppers, mixers, mills, and conveyors, should be examined in detail to increase efficiency and optimize the specification of machine design. The Distinct Element Method DEM is becoming one of the effective numerical tools to simulate dynamic behavior of granular materials. Shimizu and Cundall conducted a screw conveyor simulation using three-dimensional DEM to examine particle movement inside the conveyor, and compared the results with existing theory and empirical equations. In this paper, we focus on the effect of particle shape, in which non-spherical (i.e., two-pair) particles are used instead of spherical particles in the simulation of a horizontal screw. The results are compared to the findings with spherical particles.