Mathematical Physics - Volume II - Numerical Methods

6.5 Discrete Element Methods

259

applying Newton’s second law of motion (6.1) to each discrete element using a specified load-displacement rule to all contacts to determine new element positions. During the global conglomerate (sample) evolution, new contacts may appear (between particles that have not previously touched) and existing contacts may disappear. Therefore, the global stiffness matrix of the whole group of discrete elements must be constantly updated, from cycle to cycle. For non-cohesive (loose) materials and particle systems, there is another reason why it is necessary to update the global stiffness matrix: namely, the interactions among discrete elements, both in the normal and shear directions, are not necessarily linear, which means that the stiffnesses k n and k t , which define these contacts (e.g., Figure 6.19), must be recalculated continuously [103]-[105]. Unlike computational methods of continuum mechanics, such as FEM, in DEM the role of primary variables is played by forces and displacements. Accordingly, it is necessary to develop and apply methods for describing the continuum parameters (stresses and strains) based on these forces and displacements (i.e., mesoscopic parameters of the state of the individual elements that make up the assembly). This process is called homogenization . The starting point in this procedure is to define a representative volume element , which serves as an averaging volume for calculating mean values of macroscopic parameters (e.g., [106], [107]). Thus, homogenization and constitutive modeling techniques make it possible to take into account the micro-/meso-structure within the application of continuum methods, but only indirectly in terms of the "mean field". The DEM model parameters are typically adjusted using experimentally observed behavior and a large number of such parameters are necessary to reproduce complex phenomena. The homogenization is, very often, a demanding job because the parameters that control geometric properties and constitutive behavior do not always have a clear physical meaning and can also show complex interdependencies [10]. Despite the open questions highlighted throughout this chapter, DEM modeling (especially the mechanical behavior of geomaterials) is on the rise not only in research but also in geotechnical engineering (e.g., [108]). The main cause of this DEM popularity is the ability of natural reproduction of localization (a phenomenon ubiquitous in quasi-brittle materials with random texture) that is difficult to capture objectively by computational methods of continuum mechanics based on network discretization (e.g., FEM). 6.5.2 Contact Algorithms It is obvious from the DEM basics outlined above that contact algorithms are the essential DEM feature. Most often, realistic and detailed modeling of particle contacts is not only too complex but also unnecessary. Therefore, in the following considerations, the force of the interaction of the particles i and j is related to their overlap illustrated in Figure 6.19a. The Hertz theory of contact mechanics [109] defines the basic law of elastic contact of two spheres of radius R i and R j by a nonlinear relation f n i j ela = 4 3 1 − ν 2 i E i + 1 − ν 2 j E j ! − 1 p R i j ( u n i j ) 3 / 2 (6.64)

between the normal contact force, f n

i j , and the maximum overlap, u i j . In Equation (6.64),

( E , ν ) are pairs of elastic constants of two materials in contact and R − 1 i j = R − 1 i + R − 1 j the