Mathematical Physics - Volume II - Numerical Methods

6.5 Discrete Element Methods

267

contacts (Figure 6.20). They can act on each other exclusively through contacts which are, having in mind their simple shape, always uniquely defined. The grain overlaps are assumed to be small compared to their size to ensure that the contacts occur “at the point”. The set of meso-properties consists of the parameters of stiffness and strength of both elements and bonds. The bonds are of finite stiffness and limited strength. The force-displacement dependence in each contact relates the relative displacements of the discrete elements with the forces and moments acting on each of them. Dynamic response of the system is calculated using a finite difference algorithm comprehensively presented by Cundall and Hart [102]. The DEM simulation technique is based on the assumption that the time step is small enough so that, in one calculation cycle, the perturbations cannot extend beyond the first neighbors of each grain. Potyondy and Cundall [128] discussed in detail the advantages of an explicit numerical scheme. Since DEM is a fully dynamic formulation, attenuation can be introduced as needed to dissipate kinetic energy. This damping mimics microscopic dissipative processes in real materials, such as the internal friction or wave scattering. Figure 6.20 illustrates the way in which the bonded-particle model simulates the mechanical behavior of a group of circular grains connected by parallel bonds. The total force and moment acting in each contact consist of the contact force f i j , which is the result of the overlap of the particles (Figure 6.20a) and represents the grain behavior (Equation (6.67) with or without damping), and the force and moment, ˆ f i j and ˆ M i j which are transmitted by the parallel bond and represent the cement behavior (Figure 6.20b). These quantities contribute to the resultant force and moment acting on both circular elements (by virtue of Newton’s third law of motion) involved in the contact and represent the input data for computational integration (Newton’s second law for a dynamic system (6.1)) using an explicit finite difference scheme to obtain grain trajectories. The constitutive rule of the contacting grains is described by the same non-cohesive interaction with friction (Chapter 6.5.3 for the case without contact damping) defined with the normal and shear stiffness, k n and k t , and the friction coefficient, µ f . This contact is established as soon as the two grains overlap. The contact stiffnesses of the bond thus established (in the directions normal to, and in, the contact plane, designated respectively by superscripts n and t ) are determined by the serial connection (6.40). The overlap, although physically impermissible, mimics, in a sense, the local deformation of the grains (especially when the contact surfaces are not smooth but rough). The contact force vector of each bond can be decomposed into a normal and a shear component as already shown by (6.67) 1 . The contact behavior of the circular particle is already discussed in Chapter 6.5.3: if u n i j ≤ 0 there is a gap (note the sign convention), and the normal and shear forces are equal to zero by definition; if u n i j > 0, there is an overlap, and sliding is defined using the Mohr-Coulomb type limit (6.65). In doing so, in contrast to the normal force (which is at any time proportional to the size of the overlap with secant stiffness, k n , as the coefficient of proportionality), the shear force is calculated in an incremental manner: after establishing contact, f t in initialized to zero; from each subsequent increase in the relative displacement of the particles in the direction of the tangent, u t , there is an increase in shear force, ∆ f t = − k t ∆ u t , where k t tangential stiffness (as opposed to secant from which it is distinguished by the subscript). Contact displacements are calculated in each calculation cycle on the basis of the contact velocity which depends on the translational and angular