Implementation of the Jäger contact model for discrete element simulations

In three‐dimensional discrete element method (DEM) simulations, the particle motions within a granular assembly can produce bewildering sequences of movements at the contacts between particle pairs. With frictional contacts, the relationship between contact movement and force is non‐linear and path‐dependent, requiring an efficient means of computing the forces and storing their histories. By cleverly applying the principles of Cattaneo, Mindlin, and Deresiewicz, J”urgen Jäger developed an efficient approach for computing the full three‐dimensional force between identical elastic spheres that have undergone difficult movement sequences (J. Jäger, New Solutions in Contact Mechanics. WIT Press: Southampton, U.K.). This paper presents a complete Jäger algorithm that can be incorporated into DEM codes and also describes three special provisions for DEM simulations: (1) a method for handling particle pairs that undergo complex tumbling and twirling motions in three‐dimensions; (2) a compact data structure for storing the loading history of the many contacts in a large assembly; and (3) an approximation of the Jäger algorithm that reduces memory demand. The algorithm addresses contact translations between elastic spheres having identical properties, but it does not resolve the tractions produced by twisting or rolling motions. A performance test demonstrates that the algorithm can be applied in a DEM code with modest increases in computation time but with more substantial increases in required storage. Copyright © 2011 John Wiley & Sons, Ltd.     

The discrete element method with deformable particles

This work presents a new original formulation of the discrete element method (DEM) with deformable cylindrical particles. Uniform stress and strain fields are assumed to be induced in the particles under the action of contact forces. Particle deformation obtained by strain integration is taken into account in the evaluation of interparticle contact forces. The deformability of a particle yields a nonlocal contact model, it leads to the formation of new contacts, it changes the distribution of contact forces in the particle assembly, and it affects the macroscopic response of the particulate material. A numerical algorithm for the deformable DEM (DDEM) has been developed and implemented in the DEM program DEMPack. The new formulation implies only small modifications of the standard DEM algorithm. The DDEM algorithm has been verified on simple examples of an unconfined uniaxial compression of a rectangular specimen discretized with regularly spaced equal bonded particles and a square specimen represented with an irregular configuration of nonuniform‐sized bonded particles. The numerical results have been verified by a comparison with equivalent finite element method results and available analytical solutions. The micro‐macro relationships for elastic parameters have been obtained. The results have proved to have enhanced the modeling capabilities of the DDEM with respect to the standard DEM.     

Time reversibility of the discrete element method

An algorithm is presented for discrete element method simulations of energy-conserving systems of frictionless, spherical particles in a reversed-time frame. This algorithm is verified, within the limits of round-off error, through implementation in the LAMMPS code. Mechanisms for energy dissipation such as interparticle friction, damping, rotational resistance, particle crushing, or bond breakage cannot be incorporated into this algorithm without causing time irreversibility. This theoretical development is applied to critical-state soil mechanics as an exemplar. It is shown that the convergence of soil samples, which differ only in terms of their initial void ratio, to the same critical state requires the presence of shear forces and frictional dissipation within the soil system.     

Numerical simulation of granular materials by an improved discrete element method

In this paper, we present an improved discrete element method based on the non‐smooth contact dynamics and the bi‐potential concept. The energy dissipated during the collisions is taken into account by means of restitution coefficients. The interaction between particles is modelled by Coulomb unilateral contact law with dry friction which is typically non‐associated: during the contact, the sliding vector is not normal to the friction cone. The main feature of our algorithm is to overcome this difficulty by means of the bi‐potential theory. It leads to an easy implement predictor–corrector scheme involving just an orthogonal projection onto the friction cone. Moreover the convergence test is based on an error estimator in constitutive law using the corner stone inequality of the bipotential. Then we present numerical simulations which show the robustness of our algorithm and the various possibilities of the software ‘MULTICOR’ developed with this approach. Copyright © 2004 John Wiley & Sons, Ltd.     

A novel discrete element method based on the distance potential for arbitrary 2D convex elements

A new 2‐dimensional discrete element method, which is able to simulate a system involving a large number of arbitrary convex elements, is proposed. In this approach, a novel distance potential function is defined using a normalized format of the penetrated distance between contact couples, while a holonomic and precise algorithm for contact interaction is established, accounting for the influence of the tangential contact force. Furthermore, the new contact detection algorithm is well suited for nonuniform blocks unlike the common no binary search method that requires uniform elements. The proposed method retains the merit of the combined finite‐discrete element method and avoids its deficiencies. Compared with the existing finite‐discrete element method, the distance potential function has a clear physical meaning, where the calculation of contact interaction avoids the influence of the element shape. Accordingly, the new method completely gets rid of the restraint of uniform element type and can be applied to arbitrary convex elements. The new method is validated with well‐known benchmark examples, and the results are in very good agreement with existing experimental measurement and analytical solutions. Finally, the proposed method is applied to simulate the Tangjiashan landslide.     

Examining the smooth and nonsmooth discrete element approaches to granular matter

The smooth and nonsmooth approaches to the discrete element method (DEM) are examined from a computational perspective. The main difference can be understood as using explicit versus implicit time integration. A formula is obtained for estimating the computational effort depending on error tolerance, system geometric shape and size, and on the dynamic state. For the nonsmooth DEM (NDEM), a regularized version mapping to the Hertz contact law is presented. This method has the conventional nonsmooth and smooth DEM as special cases depending on size of time step and value of regularization. The use of the projected Gauss‐Seidel solver for NDEM simulation is studied on a range of test systems. The following characteristics are found. First, the smooth DEM is computationally more efficient for soft materials, wide and tall systems, and with increasing flow rate. Secondly, the NDEM is more beneficial for stiff materials, shallow systems, static or slow flow, and with increasing error tolerance. Furthermore, it is found that just as pressure saturates with depth in a granular column, due to force arching, also the required number of iterations saturates and become independent of system size. This effect make the projected Gauss‐Seidel solver scale much better than previously thought. Copyright © 2014 John Wiley & Sons, Ltd.     

A combined finite‐discrete element method for simulating pharmaceutical powder tableting

The pharmaceutical powder and tableting process is simulated using a combined finite‐discrete element method and contact dynamics for irregular‐shaped particles. The particle‐scale formulation and two‐stage contact detection algorithm which has been developed for the proposed method enhances the overall calculation efficiency for particle interaction characteristics. The irregular particle shapes and random sizes are represented as a pseudo‐particle assembly having a scaled up geometry but based on the variations of real powder particles. Our simulations show that particle size, shapes and material properties have a significant influence on the behaviour of compaction and deformation. Copyright © 2005 John Wiley & Sons, Ltd.     

Linear-frictional contact model for 3D discrete element simulations of granular systems

The linear-frictional contact model is the most commonly used contact mechanism for discrete element (DEM) simulations of granular materials. Linear springs with a frictional slider are used for modeling interactions in directions normal and tangential to the contact surface. Although the model is simple in two dimensions, its implementation in 3D faces certain subtle challenges, and the particle interactions that occur within a single time step require careful modeling with a robust algorithm. The paper details a three-dimensional algorithm that accounts for the changing direction of the tangential force within a time step, the transition from elastic to slip behavior within a time step, possible contact sliding during only part of a time step, and twirling and rotation of the tangential force during a time step. Without three of these adjustments, errors are introduced in the incremental stiffness of an assembly. Without the fourth adjustment, the resulting stress tensor is not only incorrect but it is also no longer a tensor. The algorithm also computes the work increments during a time step, both elastic and dissipative.     

A two-layer mesh method for discrete element simulation of gas-particle systems with arbitrarily polyhedral mesh

Inter-phase momentum coupling for particle flows is usually achieved by means of direct numerical simulation (DNS) or point source method (PSM). DNS requires the mesh size of the continuous phase to be much smaller than the size of the smallest particle in the system, whereas PSM requires the mesh size of the continuous phase to be much larger than the particle size. However, for applications where mesh sizes are similar to the size of particles in the system, neither DNS nor PSM is suitable. In order to overcome the dependence of mesh on particle sizes associated with DNS or PSM, a two-layer mesh method (TMM) is proposed. TMM involves the use of a coarse mesh to track the movement of particle clouds and a fine mesh for the continuous phase, with mesh interpolation for information exchange between the coarse and fine mesh Numerical tests of different interpolation methods show that a conservative interpolation scheme of the second order yields the most accurate results. Numerical simulations of a fluidized bed show that there is a good agreement between predictions using TMM with a second-order interpolation scheme and the experimental results, as well as predictions obtained with PSM. Copyright © 2015 John Wiley & Sons, Ltd.     

A new algorithm for contact detection between spherical particle and triangulated mesh boundary in discrete element method simulations

In discrete element method (DEM) simulations of real scale, the spherical particles are commonly employed for increasing the computation speed, and the complex boundary models are represented by triangle meshes with controllable accuracy. A new contact detection algorithm has been developed to resolve the contacts between the spheres and the triangle mesh boundaries. The application of the barycentric coordinates makes this algorithm more efficient to identify contacts in the intersection test. As a particle probably collides with several triangles at the same time, the multiple contacts would be reported as face contacts, edge contacts, or vertex contacts. Moreover, the particle embedding in a triangle can be also contact with the edges or vertices of the next triangles. These contacts should be considered as invalid for updating contact forces in the DEM. To exclude invalid records from the multiple contacts, the algorithm gives attention to the mesh structure nearby contacts and analyzes all possible collision situations. Numerical experiments have been conducted to verify this algorithm by using the algorithm in the DEM simulation framework. The numerical results suggest that the algorithm can resolve all contacts precisely and stably when the spherical particles collide on the complex boundary circumstances. Copyright © 2013 John Wiley & Sons, Ltd.     

Analytical study of the accuracy of discrete element simulations

The numerical errors was used to verify the correctness of key results. The truncation errors, which are larger than the round‐off errors by orders of magnitude, have a superlinear relationship with both the simulation time‐step and the interparticle collision speed. This remains the case regardless of the simulation details including the chosen contact model, particle size distribution, particle density or stiffness. Hence, the total errors can usually be reduced by choosing a smaller time‐step. Increasing the polydispersity in a simulation by including smaller particles necessitates choosing a smaller time‐step to maintain simulation stability and reduces the truncation errors in most cases. The truncation errors are increased by the dissipation of energy by frictional sliding or by the inclusion of damping in the system. The number of contacts affects the accuracy, and one can deduce that because 2D simulations contain fewer interparticle contacts than the equivalent 3D simulations, they therefore have lower accrued simulation errors. Copyright © 2016 John Wiley & Sons, Ltd.     

A formulation of arbitrarily shaped surface elements for three-dimensional large deformation contact with friction

The paper introduces a general theory for the numerical simulation of large deformation contact problems. The contacting bodies under consideration may be of two- or three-dimensional shape modelled by finite elements. A contact finite element which can be applied to handle multi-body contact as well as contact with rigid bodies is developed. The element is universal in the sense that it can be used as a surface element for any known finite element model and includes friction. The frictional behaviour of the model obeys Coulomb's law of friction distinguishing between sticking and sliding contact. The algorithmic treatment is based on a penalty formulation for the normal and sticking contact. The corresponding consistent tangential stiffness matrices are derived, leading to an overall quadratic convergence behaviour for the method. This feature is demonstrated in a number of representative examples. © 1997 by John Wiley & Sons, Ltd.     

A static discrete element method with discontinuous deformation analysis

For discrete element methods (DEMs), integrating the equation of motion based on Newton's second law is an integral part of the computation. Accelerations and velocities are involved even for modeling static mechanics problems. As a consequence, the accuracy can be ruined and numerous calculation steps are required to converge. In this study, we propose a static DEM based on discontinuous deformation analysis (DDA). The force of inertia is removed to develop a set of static equilibrium equations for distinct blocks. It inherits the advantages of DDA in dealing with distinct block system such as jointed rock structures. Furthermore, the critical numerical artifact in DDA, ie, artificial springs between contact blocks, is avoided. Accurate numerical solution can be achieved in mere one calculation step. Last but not the least, since the method is formulated in the framework of mathematical programming, the implementation can be easily conducted with standard and readily available solvers. Its accuracy and efficiency are verified against a series of benchmarks found in the literature.     

A new contact detection method for arbitrary dilated polyhedra with potential function in discrete element method

Contact detection significantly affects the computational efficiency of discrete element simulations, especially for irregularly shaped elements. The dilated polyhedron is constructed by the Minkowski sum of a dilated sphere and a core convex polyhedron. One of the greatest advantages of using the dilated polyhedron in contact detection lies in its ability to be solved by calculating the nearest distance between corresponding core polyhedra. The approximate envelope function (AEF) of a dilated polyhedron is formed by the weighted summation of the second-order dilated function of the polyhedral and spherical functions. The AEF can be used to represent the element in the optimization model for the contact center. Geometric calculations are then employed for the contact points on the core polyhedron, whereupon the contact detection is solved. The accuracy and stability of the proposed method by a 3-D Voronoi tessellation are validated using analytical solutions and previously published simulation results. The efficiency tests show that the speedup of the CPU-based multithread algorithm can reach 14 on a desktop. The direct shear test of the Voronoi shaped ballast is analyzed by this method. The shear stress under different vertical pressure is compared with previously published experimental and simulated results.     

An accurate and robust contact detection algorithm for particle‐solid interaction in combined finite‐discrete element analysis

When applying the combined finite‐discrete element method for analysis of dynamic problems, contact is often encountered between the finite elements and discrete elements, and thus an effective contact treatment is essential. In this paper, an accurate and robust contact detection algorithm is proposed to resolve contact problems between spherical particles, which represent rigid discrete elements, and convex quadrilateral mesh facets, which represent finite element boundaries of structural components. Different contact scenarios between particles and mesh facets, or edges, or vertices have been taken into account. For each potential contact pair, the contact search is performed in an hierarchical way starting from mesh facets, possibly going to edges and even further to vertices. The invalid contact pairs can be removed by means of two reasonable priorities defined in terms of geometric primitives and facet identifications. This hierarchical contact searching scheme is effective, and its implementation is straightforward. Numerical examples demonstrated the accuracy and robustness of the proposed algorithm. Copyright © 2015 John Wiley & Sons, Ltd.     

NSCD discrete element method for modelling masonry structures

The aim of this paper is to present a discrete element approach for the study of stonework. In the present work, a masonry structure is considered as a collection of rigid or deformable blocks, interacting together by contact laws. In this paper, we use the non‐smooth contact dynamics (NSCD) resolution method mainly used for the modelling of granular media. In the considered masonry structures we define, on an elementary cell, average local strain and stress tensors. These definitions are valid under dynamical loading of the structure, taking into account rotations. We present their use on academic and on real masonry structures. Copyright © 2005 John Wiley & Sons, Ltd.     

Finite strain,finite rotation quadratic tetrahedral element for the combined finite–discrete element method

In the past, the combined finite–discrete element was mostly based on linear tetrahedral finite elements. Locking problems associated with this element can seriously degrade the accuracy of their simulations. In this work an efficient ten‐noded quadratic element is developed in a format suitable for the combined finite–discrete element method (FEMDEM). The so‐called F‐bar approach is used to relax volumetric locking and an explicit finite element analysis is employed. A thorough validation of the numerical method is presented including five static and four dynamic examples with different loading, boundary conditions, and materials. The advantages of the new higher‐order tetrahedral element are illustrated when brought together with contact detection and contact interaction capability within a new fully 3D FEMDEM formulation. An application comparing stresses generated within two drop experiments involving different unit specimens called Vcross and VRcross is shown. The Vcross and VRcross units of ～3.5 × 10⁴kg show very different stress generation implying different survivability upon collision with a deformable floor. The test case shows the FEMDEM method has the capability to tackle the dynamics of complex‐shaped geometries and massive multi‐body granular systems typical of concrete armour and rock armour layers. Copyright © 2009 John Wiley & Sons, Ltd.     

Penalty function method for combined finite–discrete element systems comprising large number of separate bodies

Large‐scale discrete element simulations, the combined finite–discrete element method, DDA as well as a whole range of related methods, involve contact of a large number of separate bodies. In the context of the combined finite–discrete element method, each of these bodies is represented by a single discrete element which is then discretized into finite elements. The combined finite–discrete element method thus also involves algorithms dealing with fracture and fragmentation of individual discrete elements which result in ever changing topology and size of the problem. All these require complex algorithmic procedures and significant computational resources, especially in terms of CPU time. In this context, it is also necessary to have an efficient and robust algorithm for handling mechanical contact. In this work, a contact algorithm based on the penalty function method and incorporating contact kinematics preserving energy balance, is proposed and implemented into the combined finite element code. Copyright © 2000 John Wiley & Sons, Ltd.     

Validation of microparameters in discrete element modeling of sea ice failure process

A GPU-based discrete element method (DEM) with bonded particles is investigated to simulate the mechanical properties of sea ice in uniaxial compressive and three-point bending tests. Both the uniaxial compressive strength and flexural strength of sea ice are related to the microparameters in DEM simulation including particle size, sample size, bonding strength, and interparticle friction coefficient. These parameters are analyzed to build the relationship between the material macrostrengths of sea ice and the microparameters of the numerical model in DEM simulations. Based on this relationship, the reasonable microparameters can be calculated by given macrostrengths in the applications of simulating the failure processes of sea ice. In this simulation, both uniaxial compressive strength and flexural strength of ice increase with the increasing ratio of sample size and particle size. The interparticle friction coefficient is directly related to the compressive strength but has little effect on the flexural strength. In addition, numerical simulations are compared with experimental data to show the performance of the proposed model, and a satisfactory agreement is achieved. Therefore, this microparameter validation approach based on macrostrengths can be applied to simulate the complicated failure process of sea ice interacting with offshore platform structures.