Background — the 3D Distinct Element Method
3DEC is a numerical program based on the distinct element method. This method falls within the general classification of discontinuum analysis techniques. A discontinuous medium is distinguished from a continuous medium by the existence of interfaces or contacts between the discrete bodies that make up the system. Discontinuum methods can be categorized both by the way they represent contacts and by the way they represent the discrete bodies in the numerical formulation.
Aspects of Modeling a Discontinuous System
A numerical model must represent two types of mechanical behavior in a discontinuous system: (1) behavior of the discontinuities; and (2) behavior of the solid material. First, the model must recognize the existence of contacts or interfaces between the discrete bodies that make up the system. Numerical methods are divided into two groups by the way they treat behavior in the normal direction of motion at contacts. In the first group (using a soft-contact approach), a finite normal stiffness is taken to represent the measurable stiffness that exists at a contact or joint. In the second group (using a hard-contact approach), interpenetration is regarded as nonphysical, and algorithms are used to prevent any interpenetration of the two bodies that form a contact.
The choice-of-contact assumption should be made on the basis of physics rather than numerical convenience or mathematical elegance. Depending on the circumstances involved, it is possible for the same physical system to exhibit different behavior. For example, an assembly of spheres is best represented with rigid contacts when the friction coefficient is zero and the stress level is very small (see Papadopoulos 1986). However, if wave propagation is modeled through the same assembly at higher stress and friction, the contact stiffness must be taken into account in order to obtain the correct wave speed.
The preceding comments relate to the magnitude of the contact force. In addition, the contact location must be identified in the model. For point contacts (or contact almost at a point), the location of the resultant force vector clearly is at the point of contact, but where contact conditions exist over a finite surface area on both bodies, the force location is not so obvious. One assumption might be that the resultant force acts at the centroid of the interpenetration volume. Cundall (1988) suggests that the location should be regarded as an independent constitutive property, depending on the relative rotation of the two surfaces in contact. Even if a computer program can relate force location to geometric variables, there are, at present, very little data from physical tests to substantiate any physical assumption.
The second type of mechanical behavior that the model must represent is the behavior of the solid material that constitutes the particles or blocks in the discontinuous system. There are two main divisions in this representation: the material may be assumed to be rigid or deformable. The assumption of material rigidity is a good one when most of the deformation in a physical system is accounted for by movement on discontinuities. This condition applies, for example, in an unconfined assembly of rock blocks at a low stress level, such as a shallow slope in well-jointed rock. The movements consist mainly of sliding and rotation of blocks, and of opening and interlocking of interfaces.
If the deformation of the solid material cannot be neglected, two main methods can be used to include deformability. In the direct method of introducing deformability, the body is divided into internal elements or boundary elements in order to increase the number of degrees-of-freedom. The possible complexity of deformation depends on the number of elements into which the body is divided. For example, 3DEC automatically discretizes any block into tetrahedral, constant-strain zones (see Deformable Block Motion). In the elastic case, the formulation of these zones is identical to that of constant-strain finite elements. The zones may also follow an arbitrary, nonlinear constitutive law. A disadvantage of the method is that a body of complex shape must necessarily be divided into many zones, even if only a simple deformation pattern is required.
A complex deformation pattern may also be achieved in a body by the superposition of several mode shapes for the whole body. For example, Williams and Mustoe (1987) rewrite the matrix equation of motion for an element in terms of a set of orthogonal modes that may or may not be eigenmodes. In order to obtain the required complexity of deformation pattern, any number of these modes may be added. The approach is very efficient for bodies of complicated shape that deform in a simple manner, because only a few low modes need to be taken. However, it is not easy to incorporate material nonlinearity, because of the need for superposition.
A somewhat similar scheme was devised by Shi (1989) in his “discontinuous deformation analysis” (DDA). This method uses series approximations to supply an increasingly complex set of strain patterns that are superimposed for each block. However, the use of direct strain modes may be inconsistent (Williams and Mustoe 1987); the comment also applies to the “simply deformable” element of Cundall et al. (1978).
Computer Programs for Modeling Discontinuous Systems
Many computer programs based upon a continuum mechanics formulation (e.g., finite element and Lagrangian finite-difference programs) can simulate the variability in material types and nonlinear constitutive behavior typically associated with a rock mass, but the representation of discontinuities requires a discontinuum-based formulation. There are several finite element, boundary element and finite difference programs available. These programs have interface elements or “slide lines” that enable them to model a discontinuous material to some extent. However, their formulation is usually restricted in one or more of the following ways: first, the logic may break down when many intersecting interfaces are used; second, there may not be an automatic scheme for recognizing new contacts; and, third, the formulation may be limited to small displacements and/or rotation. For these reasons, continuum codes with interface elements are restrictive in their applicability for analysis of underground excavations in jointed rock.
A class of computer programs, collectively described as discrete element codes, provides the capability to represent the motion of multiple, intersecting discontinuities explicitly. Cundall and Hart (1992) provide the following definition of a discrete element method: the name “discrete element” applies to a computer program only if it:
- allows finite displacements and rotations of discrete bodies, including complete detachment; and
- recognizes new contacts automatically as the calculation progresses.
A discrete element code typically will embody an efficient algorithm for detecting and classifying contacts. It will maintain a data structure and memory allocation scheme that can handle many hundreds or thousands of discontinuities.
Cundall and Hart (1992) identify the following four main classes of codes that conform to the definition of a discrete element method:
- Distinct element programs use an explicit time-marching scheme to solve the equations of motion directly. Bodies may be rigid or deformable (by subdivision into elements); contacts are deformable. “Static relaxation” is a variation. Representative codes are TRUBAL (Cundall and Strack 1979a), UDEC (Cundall 1980; Cundall and Hart 1985; Itasca 2011), 3DEC (Cundall 1988; Hart et al. 1988), DIBS (Walton 1980), 3DSHEAR (Walton et al. 1988) and PFC (Itasca 2008).
- Modal methods are similar to the distinct element method in the case of rigid blocks, but for deformable bodies, modal superposition is used (e.g., Williams and Mustoe 1987). This method appears to be better-suited for loosely packed discontinua; in dynamic simulation of dense packings, eigenmodes are apparently not revised to account for additional contact constraints. A representative code is CICE (Hocking et al. 1985).
- Discontinuous deformation analysis assumes contacts and bodies may be rigid or deformable, depending on the version used. The condition of no penetration is achieved by an iterative scheme; the deformability comes from superposition of strain modes. The relevant computer program is DDA (Shi 1989).
- Momentum-exchange methods assume both the contacts and bodies to be rigid: momentum is exchanged between two contacting bodies during an instantaneous collision. Friction sliding can be represented (for example, see Hahn 1988).
Another class of codes, defined as limit equilibrium methods, can also model multiple intersecting discontinuities, but does not satisfy the requirements for a discrete element code. These codes use vector analysis to establish whether it is kinematically possible for any block in a blocky system to move and become detached from the system. This approach does not examine subsequent behavior of the system of blocks or redistribution of loads. All blocks are assumed rigid. The “key-block” theory by Goodman and Shi (1985) and the vector stability analysis approach by Warburton (1981) are examples of this method.
Cundall and Hart (1992) summarize the attributes of the various discrete element and limit equilibrium methods (figure below). The class of finite element or finite difference methods with slide lines is not included because of the great variations between programs. There are some programs in this class that exhibit most of the capabilities shown below, but they do not have both automatic contact detection and general interaction logic, including finite rotations and interlocking of blocks.
Figure 1: Attributes of the four classes of the discrete element method and the limit equilibrium method (Cundall and Hart 1992).
History of the Distinct Element Method
The formulation and development of the distinct element method has progressed for a period of over 40 years, beginning with the initial presentation by Cundall (1971). Figure 2 shows a chronological chart of the development of the method and relevant papers by Dr. Cundall and his associates.
The distinct element method was originally created as a two-dimensional representation of a jointed rock mass, but the method has also been extended to applications in particle flow research (see Walton 1980), studies on microscopic mechanisms in granular material (see Cundall and Strack 1983), and crack development in rocks and concrete (see Plesha and Aifantis (1983) and Lorig and Cundall (1987)). Distinct element models of jointed-rock problems have been made by many investigators (e.g., Bardet and Scott (1985), Butkovich et al. (1988), Cundall (1974), Hart et al. (1990), Heuzé et al. (1990)). The two-dimensional program, UDEC (Cundall (1980), and Lemos et al. (1985)), was first developed in 1980 to combine, into one code, formulations to represent both rigid and deformable bodies (blocks) separated by discontinuities. This code can perform either static or dynamic analyses.
In 1983, Dr. Cundall began work on the development of a three-dimensional version of the method. This work is embodied in 3DEC (Cundall 1988; Hart et al. 1988).
The most recent distinct element developments are the two-dimensional and three-dimensional particle flow codes, PFC2D and PFC3D. These codes can be applied to simulate both granular materials (such as sand) and bonded materials (such as concrete and rock). Fracturing is simulated in PFC2D and PFC3D by progressive bond breakage under load (Potyondy et al. 1996; Itasca 2008).
[cs: need the figure below updated, please]
Figure 2: Chronology of the distinct element method.
Numerical Formulation of the Distinct Element Method in 3D
This section discusses the numerical formulation for a three-dimensional distinct element model. The major elements of this formulation are the scheme for contact detection and representation in three-dimensions, and the mechanical calculations for motion and interaction in three-dimensions. This discussion is taken from the publications by Cundall 1988 and Hart et al. 1988.
A Scheme to Detect and Represent Contacts in 3D
The distinct element method has advanced to a stage where the complex mechanical interactions of a discontinuous system can be modeled in three dimensions. An important component is the formulation of a robust and rapid technique to detect and categorize contacts between three-dimensional particles. The technique, described in this section, can detect the contact between blocks of any arbitrary shape (convex or concave) and represent the geometrical and physical characteristics prescribed for the contact (e.g., three-dimensional rock joint behavior). The method utilizes an efficient data structure which permits the rapid calculation of systems involving several hundred particles on a personal computer.
The distinct element method is a way to simulate the mechanical response of systems composed of discrete blocks or particles (Cundall and Strack 1979b). Particle shapes are arbitrary: any particle may interact with any other particle, and there are no limits placed on particle displacements or rotations. In view of this generality, a robust and rapid method must be found to identify pairs of particles that are touching, and to represent their geometric and physical characteristics (e.g., whether faces, edges or vertices are involved, and what the direction of potential sliding might be).
This section describes a way to perform this task rapidly for a three-dimensional system composed of many blocks. In general, the blocks may be convex or concave, with faces that consist of arbitrary, plane polygons.
Throughout this section, reference will be made to various elements of the data structure. It is important to have a data structure that allows relevant data to be retrieved rapidly when needed, particularly in view of the explicit nature of the mechanical calculations, which often entail many thousands of passes through the main cycle. Note that Key (1986) describes a data structure for 3D sliding interfaces, but in his scheme the potential interactions must be identified in advance by the user.
The Data Structure
All data are stored in a single main array that can hold real numbers or integers, mixed in whatever way is needed. In 3DEC, real numbers and integers are stored as 64-bit words. Each physical entity (such as a block, a face or a contact) is represented by a “data element,” which is a contiguous group of two or more words. Data elements are allocated dynamically from the main array, as required, and linked to the data structure by pointers. The various types of linked lists and data elements are described in the following two sections. Data elements that are no longer needed (e.g., data from a contact that has broken) are collected together in a heap. If new data elements are needed, this heap is checked before allocating fresh memory. In a program such as 3DEC, elements come in only a small range of possible sizes. Thus, a general “garbage collector” is not required, because new data elements can frequently be allocated from discarded elements of the same length. It should be noted that linked-list schemes take very little computer time to maintain: it only requires two or three integer assignments to delete or add an item to a linked list. No reordering is necessary.
Every element in the data structure has an address in the main array: an “address” is the index in the main array of the first word of the data element. For instance, rock blocks are not referred to by sequential numbers, or by user-given numbers, but by addresses in the main array. Normally, users do not need to know about addresses, because blocks, contacts, etc., are identified by their coordinates.
A guiding principle in designing the data structure has been to reduce computer time at the expense of using more memory to store data and pointers to data. Memory is rapidly becoming plentiful and cheap, while the processing speed of computers is not increasing at the same rate. Linked lists are used extensively in 3DEC. However, a linked-list data structure is not well-suited to supercomputers that derive their speed from vector processing or from pipelining, because data are not organized sequentially in memory. In contrast, the program is placed to take advantage of a machine consisting of parallel processors connected in a 3D cubic structure. The program’s data space is partitioned naturally by the cell logic described below. Each processor could take charge of one or more cells and the particles contained in them. As particles move, they are remapped to adjacent cells — and to adjacent processors, if appropriate. Interactions (contacts) that span processor boundaries are handled by the fast data buses that interconnect processors. Each processor is assumed to have enough local memory to contain all the data associated with the 3DEC cells for which it is responsible. The explicit calculations of 3DEC are well-suited for parallel processing because, conceptually, they are already done in parallel within each step. This is not true for implicit methods, in which all elements interact numerically at each step.
Representation of Polyhedra — There are two types of polyhedral block that can be modeled by 3DEC: rigid blocks, which have six degrees-of-freedom (three translational and three rotational); and deformable blocks, which are subdivided internally into tetrahedra that have three translational degrees of freedom at each vertex (or node). Rigid blocks have plane faces that are polygons of any number of sides. Figure 6 illustrates the data structure for a rigid block. Note that each element in the data structure, although drawn separately, is embedded in the main data array, and connected by pointers, as shown. Blocks are accessed via a global pointer that gives entry to a list of all blocks in arbitrary order. The data element for each block contains a pointer that gives access to lists of vertices and faces. Each face element points to a circular list that contains addresses of the vertices that make up the face, arranged in order. It is possible, then, to access vertex data in two ways: first, the vertices can be scanned directly (for example, to update their velocities and coordinates during block motion); and second, the vertices can be accessed via each face (which is useful during contact detection). The data structure for fully deformable blocks is similar, but each original polygonal face is discretized into triangular sub-faces, in accordance with the internal discretization into tetrahedra. Each sub-face, and its associated data structure, is exactly like a regular face. However, one word in the block element points to a list of all internal tetrahedral zones. The Deformable Block Motion section provides more details on deformable blocks.
Figure 3: Simplified data structure associated with a polyhedral rigid block.
As far as the user is concerned, 3DEC accepts blocks that are concave, or even hollow or multiply connected. However, there are so many advantages to convex blocks that, within the program, concave blocks are decomposed into two or more convex blocks: one is termed a “master block”; the others are “slave blocks.” In all the logic described here, slave blocks are treated in exactly the same way as master blocks in order to take advantage of convexity. However, in the mechanical calculations, the whole block (master and slaves) is treated as one. Thus, a common center of gravity, a common mass, etc., are determined. In what follows, convexity will be invoked to justify several procedures, but nothing limits the program’s ability to treat blocks of arbitrary shape.
Representation of Contacts — For each pair of convex, rigid blocks that touches (or is separated by a small enough gap), a single data element is assigned. This element corresponds to the physical contact between the two blocks, and contains relevant data such as friction, shear force, and so on. Each contact is discretized into sub-contacts, where interaction forces are applied. Sub-contacts are created for both rigid and deformable blocks. This logic is described in the Interaction between Blocks section. The single contact element contains a pointer to a series of sub-contacts.
The major items in each contact element are shown in Figure 4. Each contact is linked globally to all other contacts, as well as locally to the pair of blocks that make up the contact. The form of the data structure that embodies this linkage is shown in Figure 5, for an example system consisting of four blocks. Each contact can be accessed in several ways, depending on the need. During the main calculation cycle, when all forces are updated, it is convenient to scan through all contacts one at a time. This is done by using the linked list that is attached to the global pointer. Once a contact has been accessed, its constituent blocks are identified from the contact’s block pointers. During contact detection and reassignment, it is convenient to know what contacts with a given block already exist. Local lists, which originate with each block, thread through the block’s contacts. For example, if the pointer from block C is followed (Figure 5), the three contacts of block C are discovered.
Figure 4: Main items in contact data element.
Figure 5: Local and global contact links for a four-block example.
Identification of Neighbors
Before the relative geometry of a pair of blocks can be investigated by the computer program, candidate pairs must be identified. It is prohibitive, in computer time, to check all possible pairs, as the search time increases quadratically with the number of blocks. In two dimensions, it is possible to set up a “canonic” data structure that represents the voids between blocks automatically (Cundall 1980). It is then a simple matter to scan the local voids surrounding a block, in order to obtain a list of all possible contacting blocks. This scheme has a search time that increases linearly with the number of blocks in a system, but unfortunately, the data structure does not translate into three dimensions. A less elegant, but perhaps more robust, scheme was used in the programs BALL and TRUBAL, which model disks and spheres, respectively (Cundall and Strack 1979a). This is the method adopted here for the identification of neighbors.
Cell Mapping and Searching — The space containing the system of blocks is divided into rectangular 3D cells. Each block is mapped into the cell or cells that its “envelope space” occupies. A block’s envelope space is defined as the smallest three-dimensional box with sides parallel to the coordinate axes that can contain the block. Each cell stores, in linked-list form, the addresses of all blocks that map into it. Figure 6 illustrates the mapping logic for a two-dimensional space (as it is difficult to illustrate the concept in three dimensions). Once all blocks have been mapped into the cell space, it is an easy matter to identify the neighbors to a given block: the cells that correspond to its envelope space contain entries for all blocks that are near. Normally, this “search space” is increased in all directions by a tolerance, so that all blocks within the given tolerance are found. Note that the computer time necessary to perform the map and search functions for each block depends on the size and shape of the block, but not on the number of blocks in the system. The overall computer time for neighbor detection is consequently directly proportional to the number of blocks, provided that cell volume is proportional to average block volume.
Figure 6: Examples of block mapping to cell space, in two dimensions.
It is difficult to provide a formula for optimum cell size because of the variety of block shapes that may be encountered. In the limit, if only one cell is used, all blocks will map into it, and the search time will be quadratic. As the density of cells increases, the number of non-neighboring blocks retrieved for a given block will decrease. At a certain point, there is no advantage to increasing the density of cells, because all the blocks retrieved will be neighbors. However, by further increasing the cell density, the time associated with mapping and searching increases. The optimum cell density must therefore be of the order of one cell per block, in order to reduce both sources of wasted time.
Scheme for Triggering Neighborhood Searchers — As a block moves during the course of the simulation, it is remapped and tested for contact with new neighbors. This process is triggered by the accumulated movement of the block — a variable \(u_{acc}\), set to zero after each remap, is updated at every timestep, as follows:
where \(du\) is the incremental displacement of a vertex, and the \(\max\{ \}\) function is taken over all vertices of the block.
When \(u_{acc}\) exceeds CTOL unsure what font CTOL should be (CTOL is a preset tolerance that can be viewed with the command block list information and changed with the command block contact tolerance), remapping and contact testing are activated. The contact testing is done for a search volume that is 2\(\cdot\)CTOL larger in all dimensions than the block envelope. In this way, maximum movement of the block, and any potential neighbor, is allowed. If any block attempts to move outside the cell space (i.e., the total volume covered by cells), the cell space is redefined to be 10% larger in the affected dimension. In this case, a complete remap of all blocks occurs.
The value of CTOL is also used to determine whether a contact is created or deleted. If two blocks are found to be separated by a gap that is equal to or less than CTOL, a contact is created. Conversely, if an existing contact acquires a separation that is greater than CTOL, the contact is deleted. The logic described above ensures that the data structure for all potential contacts is in place before physical contact takes place. It also ensures that contact searching is only done for moving blocks; there is no time wasted on relatively inactive blocks.
Contact Detection
Requirements of the Scheme — After two blocks have been recognized as neighbors, they are then tested for contact: if they are not in contact, the maximum gap between them must be determined so that block-pairs separated by more than a certain tolerance may be ignored. For block-pairs separated by less than this tolerance, but not touching, a “contact” is still formed. Though the contact carries no load, it is tracked at every step in the mechanical calculation. In this way, interaction forces start to act as soon as the blocks touch. (Note that contact detection, a lengthy process, is not done at every mechanical step.) The contact-detection logic must also supply a unit normal vector, which defines the plane along which sliding can occur. This unit normal should be well-behaved (i.e., it should change direction in a continuous fashion) as the two blocks move relative to one another. The logic should be able to handle, in a reasonable way, certain extreme cases, such as that illustrated here.
Figure 7: Extreme case in which the determination of contact normal is difficult.
Finally, the contact-detection logic must classify the type of the contact rapidly (e.g., face-to-edge or vertex-to-face). This information is needed in order to select the most appropriate physical law to apply at each contact. In summary, the contact-detection logic must supply, with as little delay as possible, the contact type (if touching), the maximum gap (if not touching), and the unit normal vector. A direct approach is described first. The difficulties with this approach are pointed out, and a better scheme is described.
Direct Tests for Contact between Two Blocks — The simplest approach is to test all possibilities for interaction. In three dimensions, there are many ways for blocks to touch one another (e.g., each vertex of the first block may be tested for contact with each vertex, edge, and face of the second block, and so on). If the first block (block A) has \(v_A\) vertices, \(e_A\) edges, and \(f_A\) faces, and the second block (block B) has \(v_B\) vertices, \(e_B\) edges, and \(f_B\) faces, the number of distinct contact combinations is
As an example, two cubes give rise to 676 contact possibilities. In practice, not so many tests are needed, because some types of contact are encompassed by other types, as limiting cases. It appears that only vertex-to-face and edge-to-edge contacts need to be checked, as the other types of contact may be recovered in the following ways:
- vertex-to-vertex is detected when three or more vertex-to-face contacts exist at the same location;
- vertex-to-edge is recognized when two vertex-to-face contacts coincide;
- edge-to-face occurs when two edge-to-edge contacts exist between two blocks; and
- face-to-face is recognized by three or more edge-to-edge contacts, or three or more vertex-to-face contacts.
Even when this reduced set of tests is done, the number of combinations is
For the two-cube example, this number is 240.
Two observations are relevant to what will follow. First, the number of tests depends quadratically on the average number of block edges (or vertices or faces). Second, the tests are not always simple: e.g., in a vertex-to-face test, it is not sufficient to just check whether the vertex lies above or below the face; it is also necessary to show that the vertex lies within the projected perimeter of the polygon that makes up the face. This is not easy to do quickly.
Considering the requirements noted above, the contact type is determined in the course of the detection calculation, as each type is checked in turn. The determination of the unit normal is easy in some cases (e.g., in all cases involving a face), but difficult in others (in particular, for edge-to-edge, edge-to-vertex, and vertex-to-vertex contacts, when undefined). Furthermore, there is no guarantee that the contact normal will evolve in a smooth way when there is a jump from one contact type to another. Using this scheme of direct testing, the determination of maximum gap between two arbitrary, non-touching blocks is not a trivial task. In response to these difficulties, the following new scheme was conceived.
The Idea of a Common Plane — The difficulties noted in the previous section arise from the need to test one arbitrary polyhedron directly with another. Many of the difficulties would vanish if the problem could be split into the following two parts: (1) determining a “common-plane” that, in some sense, bisects the space between the two blocks; and (2) testing each block separately for contact with the common-plane. The “common-plane” is analogous to a metal plate that is held loosely between the two blocks (see Figure 8). If the blocks are held tightly and brought together slowly, the plate will be deflected by the blocks and will become trapped at some particular angle when the blocks finally come into contact.
Figure 8: Visualization for positioning of common-plane in response to block geometry.
Whatever the shape and orientation of the blocks (provided they are convex), the plate will take up a position that defines the sliding plane for the two blocks. To carry the analogy a bit further, imagine that the plate is now repelled by the blocks even when they do not touch. As the blocks are brought together, the plate will take up a position midway between them, at a maximum distance from both. Then we can easily find the gap between the blocks, simply by adding the block-to-plate distances. In fact, there are many things that would become easier if we could somehow contrive a numerical equivalent for the metal plate (the common-plane referred to above). Suppose for the moment that we do have a way to do this. The task of testing for contact is simplified and speeded up in the following ways. (The term common-plane is abbreviated below as “c-p.”)
Only simple vertex-to-plane tests need to be carried out (using dot products). Because the blocks (or sub-blocks) are convex (see above), face and edge contacts are recognized simply by counting the number of vertex-to-plane overlaps for both blocks.
The number of tests depends linearly on the number of vertices (compared to the quadratic dependence noted above). Because we test the vertices of block A with the c-p, and separately test the vertices of block B with the c-p, the number of tests is
(4)\[n = v_A + v_B\]For the two-cube example, the number of tests is 16, as compared with the 240 noted above.
There is no need to test whether a potential contact lies within the perimeter of a face. If both blocks touch the c-p, then they must touch each other. (If they did not touch, then the c-p would touch neither, as the c-p is defined as bisecting the space between blocks.)
The unit contact normal is equal to the c-p unit normal — no additional calculations are necessary.
Since the c-p normal is uniquely defined, the problem of discontinuous evolution of the contact normal is eliminated. The c-p normal may change rapidly (as in vertex-to-vertex contact), but it will not jump due to a change in contact type.
The determination of minimum gap between two non-touching blocks is trivial: it is simply the sum of the distances of the two blocks from the c-p.
Algorithm to Position the Common-Plane — Having established that the common-plane simplifies and speeds up contact testing, we need to provide a means to position the common-plane and to show that the overhead associated with it does not outweigh the benefits that it brings to contact testing. The algorithm for locating and moving the c-p is based on geometry alone, and is applied every 10 steps (by default), in parallel with the mechanical calculations. The algorithm is stated:
Maximize the gap between the c-p and the closest vertex.
Figure 9 shows several examples of c-p in two dimensions, which are consistent with the algorithm given above. Any rotation or shift of the c-p would reduce the gap between the c-p and the closest vertex, or leave it unchanged. For overlapping blocks, the same algorithm applies, but the words “gap” and “closest” must be used in their mathematical sense for the case of negative signs (i.e., gap means “negative overlap” and closest means “most deeply buried”). To improve readability, the algorithm may be restated for the case of overlapping blocks:
Minimize the overlap between the c-p and the vertex with the greatest overlap.
Starting conditions for two blocks which are not already recognized as being in contact must be provided to the algorithm. An initial guess is made: the c-p is placed midway between the centroids of the two blocks, with a unit normal vector pointing from one centroid to the other:
where:
\(Z_i\) = \(B_i - A_i\);
\(z^2\) = \(Z_i\)\(Z_i\);
\(n_i\) = unit c-p normal;
\(C_i\) = reference point of c-p;
\(A_i\) = position vector of block A’s centroid;
\(B_i\) = position vector of block B’s centroid; and
indices \(i\), \(j\), and \(k\) take the values 1 to 3, and denote components of a vector or tensor. (The summation convention applies for repeated indices.)
Figure 9: Examples of the common-plane between two blocks.
The algorithm then applies a translation and a rotation to the c-p in order to maximize the gap (or minimize the overlap). The function of the reference point, \(C_i\), is twofold: it is the point about which c-p rotations are applied; and it is the reaction point for normal and shear forces when the two blocks are touching.
Translation of the Common-Plane — The translation is divided into parts that are normal and tangential to the c-p. The normal translation is found by scanning each block for its nearest vertex to the c-p:
where:
\(d_A\) is the distance to the nearest vertex on \(A\) (negative for a gap);
\(d_B\) is the distance to the nearest vertex on \(B\) (positive for a gap);
\(V_i\) (\(A\)) is the position vector of a vertex on \(A\);
\(V_i\) (\(B\)) is the position vector of a vertex on \(B\);
\(\min\){ } is the minimum, taken over all vertices of \(B\); and
\(\max\){ } is the maximum taken over all vertices of \(A\).
The shift in the c-p’s reference point is then
The total gap is \(d_B - d_A\). If the blocks are touching (i.e., the total gap is negative), the reference point becomes the reaction point for contact forces and is determined in the part of the program that deals with the mechanical calculation (see below). For non-touching blocks, the reference point is moved midway between the nearest vertices:
where:
\(V_i\)(\(A_{\max}\)) is the vertex on A nearest to the c-p; and
\(V_i\)(\(B_{\min}\)) is the vertex on B nearest to the c-p.
Rotation of the Common-Plane — Whereas translation of the c-p can be done in one step, rotation must be done iteratively, because the nearest vertex on a block can change as the c-p is rotated. Since the unit normal can rotate about two independent axes, the maximization of the gap is equivalent to a hill-climbing process. Two orthogonal axes are chosen arbitrarily — both are orthogonal to the c-p unit normal vector. The unit normal is perturbed in each of these directions, both in a positive and negative sense, making four perturbations in all. If \(p_i\) and \(q_i\) are the orthogonal unit normal vectors, the four perturbations to \(n_i\) are:
where \(z^2 = 1 + k^2\), and \(k\) is the parameter that determines the size of the perturbation.
When a contact is being created for the first time, the parameter \(k\) is initially set to \(k_{\max}\), a value that corresponds to a five-degree perturbation in angle. The iteration procedure shown in Figure 10 is then used to search for a maximum value of the gap. Note that maximum gap is defined as \(\max\){\(d_B\), -\(d_A\)}, with \(d_B\), \(d_A\) given by Equations (6) and (7).
Figure 10: Iteration procedure to search for the maximum gap.
The iteration stops when the gap decreases for all four of the smallest perturbations; it follows that this final unit normal is the one that gives rise to the largest gap. The smallest perturbation, \(k_{\min}\), has a value that corresponds to an angle change of 0.01 degree. In order to prevent the iteration from terminating prematurely on a saddle point, the perturbation axes are rotated by 45 degrees on alternate cycles of the iteration. If at any stage in the iteration, the gap exceeds the tolerance set for contact formation (CTOL), the iteration process halts and the contact is deleted.
When a contact already exists, the four perturbations are initially tried with \(k\) = \(k_{\min}\). If there is no increase in the maximum gap, nothing further is done. Otherwise, the iteration given above is performed.
Translation of the c-p during the Mechanical Cycle — When forces exist on a contact, the normal translation of the c-p is still done according to Equation (8). However, the following two further incremental translations are applied to the c-p.
(1) Rigid Body Translation
The first part of the additional translation is related to the average motion of the two blocks, at the contact point:
where \(du_i\)(\(A\)) and \(du_i\)(\(B\)) are the incremental displacements of blocks A and B, respectively, at the contact point (reference point).
For example, if both blocks are moving through space at the same speed, the reference point on the c-p will also be moved at this speed. In such a case, the correction supplied by Equation (8) will be unnecessary. Note that the incremental displacements in Equation (11) include the effects of block rotation.
(2) Relative Rotation
As mentioned before, the reference point of the c-p is assumed to be the point at which the resultant contact force acts. As the upper surface of a contact plane rotates relative to the lower surface, the resultant contact force will move, since the plane will become unequally loaded. Figure 11 illustrates this effect. The relation between the movement of the reference point and the rotation of the surfaces must depend on the nature of the interface, and on the current normal stress if the normal stiffness is stress-dependent. It appears that the relation described above is a material property and cannot be derived from geometrical factors alone. In 3DEC, therefore, the translation of the reference point is taken as the relative angle change of the two blocks multiplied by a user-specified constant, \(K_T\) :
where:
\(dT_k\) (\(A\)) is the incremental rotation vector of block A;
\(dT_k\) (\(B\)) is the incremental rotation vector of block B; and
\(e_ijk\) is the permutation tensor.
\(K_T\) , as defined above, has the dimension of length, but it should probably be normalized by the length of the contact plane in the direction of movement of \(C_i\) . However, before \(K_T\) can be defined properly, more data are needed from laboratory tests.
Figure 11: Illustration of the movement of the resultant contact force due to rotation of the surfaces in contact (overlaps are exaggerated)
(3) Limit to Movement of Reference Point
Because the reference point is the point at which contact forces act, it must lie on the surface of both blocks. After applying Equations (11) and (12), the reference point is tested against each face of the two blocks that make up the contact. If it is found to lie outside any one face, it is brought back toward the face, as follows:
where:
\(n_i\) (\(f\)) is the outward unit normal of the face; and
\(d\) is the normal distance from the face to \(C_i\).
Equation (13) brings \(C_i\) completely back to the face only when \(n_i\) is orthogonal to \(n_i\) (\(f\)). But, in other cases, the repeated use of the formula (which happens automatically, as it is invoked in every calculation cycle) achieves convergence. Figure 12 illustrates this. The same effect will occur if \(C_i\) should fall outside two or more faces simultaneously, which may happen at a corner or edge.
Figure 12: Procedure to bring reference point Ci within block boundary.
Overhead Associated with Common-Plane — The number of operations necessary to establish the c-p for a pair of blocks depends linearly on the number of vertices. For the translation correction of Equation (8), the number of tests is \(v_A\) + \(v_B\). For the rotation iteration, the number of tests is \(4N\) (\(v_A\) + \(v_B\)), where \(N\) is the number of iterations. The total number of tests is therefore
It is difficult to directly compare this formula with Equations (1) and (2), which correspond to the direct-testing scheme to find contacts. For contacts already established, the c-p scheme is superior to the direct-testing scheme, because only one rotation iteration is usually needed to confirm that the c-p position is optimal. However, on initial contact formation, the number of iterations may be in the range 9 to 30 (found by experimentation). For blocks with few vertices, the c-p scheme takes more tests to establish contact conditions than the direct-testing scheme. However, when the six previously noted advantages are considered, the c-p scheme is preferred overall.
Contact Types — Contact type is important because it determines the mechanical response of the contact. For example, an edge-to-face contact will behave differently than a face-to-face contact. In rock mechanics, face-to-face contacts are thought of as “joints” in which stresses, rather than forces, are the important variables. Contacts may be classified into types by noting how many vertices of each block touch the c-p. Table 1 relates the contact type to the number of touching vertices from each block.
For a face-to-face contact, it is necessary to define an area of contact in order to use a stress-displacement law for the mechanical behavior of the interface. Because both faces that make up the contact are convex polygons, the common area of contact is also a convex, simply connected polygon. The calculation of common area is then a straightforward, but lengthy, procedure.
In 3DEC, interaction between face-to-face contacts is represented by points of the vertex-to-face or edge-to-edge type (see Interaction between Blocks).
| Number of Vertices Touching | Contact Type | |
|---|---|---|
| Block A | Block B | |
| 0 | 0 | null |
| 1 | 1 | vertex-vertex |
| 1 | 2 | vertex-edge |
| 1 | >2 | vertex-face |
| 2 | 1 | edge-vertex |
| 2 | 2 | edge-edge |
| 2 | >2 | edge-face |
| >2 | 1 | face-vertex |
| >2 | 2 | face-edge |
| >2 | >2 | face-face |
Interaction between Blocks
The scheme described above applies for both rigid and deformable blocks — only one common-plane is found for each block pair that is in contact, and only one regular data element is allocated for the contact. If a block face is in contact with the c-p, then it is automatically discretized into sub-contacts. For rigid blocks, faces are triangulated to create the sub-contacts. These sub-contacts are generally created at the vertices of the block face. For deformable blocks, the triangular faces of tetrahedral zones at the block surface contain a number of internal surface nodes, each of which has three independent degrees of freedom. In this case, a sub-contact is created for each node on the face.
The sub-contact keeps track of the interface forces between blocks, as well as other conditions such as sliding and separation. Two types of sub-contact are defined: vertex-to-face and edge-to-edge. In order to simulate face-to-face contact, each sub-contact is assigned an area allowing standard joint constitutive relations, formulated in terms of stresses and displacements, to be applied (e.g., elastic plus Coulomb friction in the shear direction — see the section on the Coulomb-Slip Joint Model). Edge-to-edge sub-contacts model both edge-to-edge contact between blocks, and face-to-face and face-to-edge contacts at the points of intersection of edges on the c-p. The interface displacement at each sub-contact is taken as the sub-contact displacement minus the displacement of the coincident point on the opposing face. The area “owned” by each sub-contact is, in general, equal to one-third of the area of the surrounding triangles, but this calculation must be adjusted when the sub-contact is close to one or more edges on the opposing block.
We have discussed only the sub-contacts on one side of an interface. If the other side of the interface is also a face, then identical conditions apply: sub-contacts are created, and relative displacements, and hence forces, are calculated. When two blocks come together, the contact logic described above is equivalent to two sets of contact springs in parallel — in this case, the forces from both sets are divided by two, so that the overall interface behavior is the average of that of both sets.
The c-p logic described previously is strictly applicable only to convex blocks with planar faces, but these conditions may be violated if large strains occur with deformable blocks. In practice, the program is used to model a rock mass, where displacements may be large, but strains are usually quite small. In these circumstances, the logic will still work. But in situations where block strains become large (e.g., >1%), the scheme may need to be modified.
At present, 3DEC does not allow the use of rigid and deformable blocks in the same problem. The logic applies for both small-displacement and large-displacement relative motion between blocks. In order to allow large displacements, the logic incorporates a procedure to automatically relocate each sub-contact, as the associated vertex crosses a face boundary in the other block. Sub-contact locations and weights are updated every 10 steps (by default). Detection of new sub-contacts and sub-contact type changes are also performed with the same periodicity. The logic also allows the user to avoid abrupt deletion of a sub-contact whose associated vertex slides out of the other blocks’ faces. The existing sub-contact forces are reallocated to ensure a smooth transition between neighboring states, as in the two-dimensional code UDEC.
Mechanical Calculations for Motion and Interaction in 3D
This section describes the formulation of the mechanical calculations of the distinct element method in three dimensions. This discussion is based upon work presented previously by Cundall and Strack (1979b) and Cundall and Hart (1985). The description of the formulation is first given for the calculation of the mechanical interaction between blocks. This scheme, described in the Sub-contact Force Update and Coulomb-Slip Joint Model sections, applies for both rigid blocks and deformable blocks.
Following this, the Rigid Block Motion section addresses the mechanical calculation for rigid block motion. This description of motion is a sufficient representation for stability studies in which the applied stress state is low compared to the intact rock strength, and in which motion is concentrated along structural features.
3DEC also accounts for block deformability and failure of intact material. In this formulation, each polyhedral block is subdivided into an internal finite difference mesh consisting of constant-strain tetrahedral elements. The solution is an explicit, large-strain one, with elastic and elasto-plastic models for the block material. The deformable block logic is analogous to the formulation for two-dimensional blocks, as given by Lemos (1987). The description of the formulation for motion of deformable blocks is given in the Deformable Block Motion section.
Calculation Cycle
3DEC is based on a dynamic (time-domain) algorithm that solves the equations of motion of the block system by an explicit finite difference method. A solution scheme based on the equations of motion is demonstrated (Cundall 1987) to be better-suited to indicate potential failure modes of discontinuum systems than schemes which disregard velocities and inertial forces (e.g., successive over-relaxation). At each timestep, the law of motion and the constitutive equations are applied. For both rigid and deformable blocks, sub-contact force-displacement relations are prescribed. The integration of the law of motion provides the new block positions, and therefore the contact-displacement increments (or velocities). The sub-contact force-displacement law is then used to obtain the new sub-contact forces, which are to be applied to the blocks in the next timestep. The cycle of mechanical calculations is illustrated below.
Figure 13: Calculation cycle.
Sub-contact Force Update
The A Scheme To Detect And Represent Contacts In 3D section described the procedure for updating the geometric parameters associated with contact between blocks. The unit normal to the c-p is taken as the contact normal, \(n_i\) (pointing from block A to block B), and is assumed to be the same for all sub-contacts.
The relative velocity across a sub-contact is obtained from the velocity associated with the subcontact, \(V^V_i\), and the velocity of the corresponding point on the opposing face, \(V^F_i\). For rigid blocks, the velocity is calculated from the equations of motion for the rigid block (see Equation (49)). For deformable blocks, the velocity of the tetrahedral-zone vertex located on the block face is calculated with Equation (57).
For rigid blocks, the contact velocity (defined as the velocity of block B relative to block A at the sub-contact location) is calculated as
where:
\(A_i\) and \(B_i\) are the position vectors of the centroids of blocks A and B;
\(\dot x^A_i\) and \(\dot x^B_i\) are translation velocity vectors of blocks A and B;
\(\omega^A_j\) and \(\omega^B_j\) are the corresponding angular velocity vectors;
\(e_ijk\) is the permutation tensor; and
indices \(i, j, k\) take the values 1 to 3, and denote components of a vector or tensor in the global coordinate system. (The summation convention applies for repeated indices.)
The velocities of the sub-contacts associated with the rigid-block contact are then interpolated from the contact velocity.
For deformable blocks, the velocity, \(V^F_i\), can be calculated by linear interpolation of the velocities of the three vertices of the face:
The weighting factors can be calculated by transforming the coordinates of the vertices \(a\), \(b\), and \(c\) into a local reference system, with one axis normal to the face plane. Denoting the local in-plane coordinates of vertex \(a\) by \(X^a\) and \(Y^a\), the weighting factor, \(W_a\), is given by
The other two factors can be obtained by circular permutation of the superscripts \(a\), \(b\), and \(c\).
The unit normal, \(n_i\), points from block A to block B. Therefore, the relative velocity is calculated as
when the vertex belongs to block A. Otherwise it is calculated as
The increment in relative displacement at the sub-contact for both rigid and deformable blocks is given by
which can be resolved into normal and shear components along the c-p. The normal displacement increment is then given by
and the shear displacement increment vector (expressed in global coordinates) by
Note that the unit normal to the c-p, \(n_i\) , is updated at every timestep. In order to account for the incremental rotation of the c-p, the vector representing the existing shear force \(F^s_i\) (in global coordinates) must be corrected as
where \(n^{old}_m\) is the old unit normal to the c-p.
The sub-contact displacement increments are used to calculate the elastic force increments. The normal force increment, taking compressive force as positive, is
and the shear force vector increment is
where \(A_c\) = area of the sub-contact. The sub-contact area is obtained by assigning to the point a region formed by 1/3 of the areas of the triangular faces containing the sub-contact and lying on the c-p. The area of the intersection of this region with the other block’s faces lying on the c-p is then calculated.
For face-to-face contact, \(A_c\) is taken as one-half of this area, in order to account for the fact that the sub-contacts are established for vertices of both blocks, therefore resulting in two sets of parallel “springs.”
The total normal force and shear force vectors are updated for the sub-contact as
and
and adjusted according to the contact constitutive relations. The basic constitutive model in 3DEC is defined as a Coulomb-slip joint model, described in the Coulomb-Slip Joint Model section.
The sub-contact force vector, which represents the action of block A on block B, is given by
For rigid blocks, the sub-contact forces are then added to the forces and moments acting on the centroids of both blocks. The force and moment sums of block A are therefore updated as
where \(c_i\) is the position vector of the sub-contact. Similarly, for block B:
For deformable blocks, on the vertex side of the sub-contact, this force is added to the other gridpoint forces. On the face side, the force is distributed among the three vertices (\(a, b, c\)), using the interpolation factors defined above:
Coulomb-Slip Joint Model
The basic joint constitutive model incorporated in 3DEC is the generalization of the Coulomb friction law. This law works in a similar fashion both for sub-contacts between rigid blocks and for sub-contacts between deformable blocks. Both shear and tensile failure are considered, and joint dilation is included.
In the elastic range, the behavior is governed by the joint normal and shear stiffnesses, \(K_n\) and \(K_s\), as described above by Equations (24) and (25).
For an intact joint (i.e., without previous slip or separation), the tensile normal force is limited to
where \(T\) is the joint tensile strength.
The maximum shear force allowed is given by
where \(c\) and \(\phi\) are the joint cohesion [stress] and friction angle.
Once the onset of failure is identified at the sub-contact, in either tension or shear, the tensile strength and cohesion are taken as zero:
This instantaneous loss of strength approximates the “displacement-weakening” behavior of a joint. The new contact forces are corrected in the following manner (note that normal compressive force is positive):
for tensile failure
for shear failure
where the shear force magnitude, \(F^s\), is given by
Dilation takes place only when the joint is at slip. The shear increment magnitude, \(\Delta U^s\), is given by
This displacement leads to a dilation of
where \(\psi\) is the dilation angle.
The normal force must be corrected to account for the effect of dilation:
Real joints display a reduction in the dilation angle as the residual friction state is approached. In 3DEC, the joints can be prevented from dilating indefinitely by prescribing a limiting shear displacement, \(U^s_{lim}\). When the magnitude of the shear displacement exceeds \(U^s_{lim}\), the dilation angle is set to zero.
Dilation is a function of the direction of shearing. Dilation increases if the shear displacement increment is in the same direction as the total shear displacement; it decreases if the shear increment is in the opposite direction.
This joint model is illustrated below for the case that joint cohesion is initially zero.
Figure 14: Coulomb slip model (for zero joint cohesion).
A more comprehensive displacement-weakening model is also available in 3DEC. This model (the continuously yielding joint model) is intended to simulate the intrinsic mechanism of progressive damage of the joint under shear. The model is described in “Continuously Yielding Joint Model In 3DEC”.
Rigid Block Motion
The equations of translational motion for a single block can be expressed as
where:
\(\ddot x_i\) = the acceleration of the block centroid;
\(\dot x_i\) = the velocity of the block centroid;
\(\alpha\) = the viscous (mass-proportional) damping constant;
\(F_i\) = sum of forces acting on the block (from block contacts and applied external forces);
\(m\) is the block mass; and
\(g_i\) is the gravity acceleration vector.
The rotational motion of an undamped rigid body is described by Euler’s equations, in which the motion is referred to the principal axes of inertia of the body:
where:
\(I_1\), \(I_2\), \(I_3\) are principal moments of inertia of the block;
\(\dot \omega_1\), \(\dot \omega_2\), \(\dot \omega_3\) are angular accelerations about the principal axes;
\(\omega_1\), \(\omega_2\), \(\omega_3\) are angular velocities about the principal axes; and
\(M_1\), \(M_2\), \(M_3\) are components of torque applied to the block referred to the principal axes.
Rigid block models are more appropriate for quasi-static analyses, and, in these cases, the rotational equations of motion can be simplified. Because velocities are small, the nonlinear term in the preceding equations can be dropped, uncoupling the equations. Also, because the inertial forces are small compared with the total forces applied to the blocks, an accurate representation of the inertia tensor is not essential. In 3DEC, therefore, only an approximate moment of inertia, \(I\), is calculated, based upon the average distance from the centroid to vertices of the block. This allows the preceding equations to be referred to the global axes.
Inserting a viscous damping term, the equations become
where the velocities, \(\omega_i\) , and the total torque, \(M_i\) , are now referred to the global axis.
A central finite-difference procedure is used to integrate the equations of motion. The following expressions describe the translational and rotational velocities at time \(t\) in terms of the values at mid-intervals:
The accelerations are calculated as
Inserting these expressions in the equations of translational and rotational motion, Equations (44) and (46) respectively, and solving for the velocities at time \([t + (\Delta t/2)]\), results in
where:
\(D_1 = 1 - (\alpha {\Delta t \over 2})\); and
\(D_2 = {1 \over {1 + \alpha {\Delta t \over 2}}}\) .
The increments of translation and rotation are given by
The position of the block centroid is updated as
The new locations of the block vertices are given by
For groups of joined blocks, the motion calculations are only performed for the master block, whose mass, moment of inertia, and centroid position are modified to represent the group of blocks. Once the motion of the master block is determined, the new position of the centroid and vertices of the slave blocks are calculated by expressions similar to (52).
The force and moment sums, \(F_i\) and \(M_i\), for all blocks are reset to zero every cycle after the block motion update is completed.
Deformable Block Motion
For many applications, the deformation of individual blocks cannot be reasonably ignored (i.e., blocks cannot be assumed to be rigid). Fully deformable blocks were developed in 3DEC to permit internal deformation of each block in the model.
Deformable blocks are internally discretized into finite-difference tetrahedral elements. The complexity of deformation of the blocks depends on the number of elements into which the blocks are divided. The use of tetrahedral elements eliminates the problem of hourglass[1] deformations that may occur with constant-strain finite-difference polyhedra. (The term “hourglassing” comes from the shape of the deformation pattern of elements within a mesh. For polyhedra with more than four nodes, combinations of nodal displacements that produce no strain and result in no opposing forces exist. The resulting effect is unopposed deformations of alternating direction.)
The vertices of the tetrahedral elements are gridpoints, and the equations of motion for each gridpoint are formulated as
where:
\(s\) is the surface enclosing the mass, \(m\), lumped at the gridpoint;
\(n_j\) is the unit normal to \(s\);
\(F_i\) is the resultant of all external forces applied to the gridpoint (from block sub-contacts or otherwise); and
\(g_i\) is the gravitational acceleration.
Gridpoint forces are obtained as a sum of three terms:
\(F_i^l\) are the external applied loads. \(F_i^c\) result from the sub-contact forces, and exist only for gridpoints along the block boundary. Forces from sub-contacts along the two faces adjacent to the gridpoint contribute to this term. Because a linear variation of displacements is assumed along any face, the effect of sub-contact forces applied along a face may be represented by statically equivalent forces applied to the face endpoints. Finally, the contribution of the internal stresses in the zones adjacent to the gridpoint is calculated as
where:
\(\sigma_{ij}\) is the zone stress tensor; and
\(n_j\) is the unit outward normal to the contour, \(C\), which follows the closed polygonal surface defined by the straight segments which bisect the zone faces converging on the gridpoint under consideration.
A net nodal force vector, \(\sum\ F_i\), is calculated at each gridpoint. This vector includes contributions from applied loads, as discussed above, and from body forces due to gravity. Gravity forces, \(F_i^{(g)}\), are computed from
where \(m_g\) is the lumped gravitational mass at the gridpoint, defined as the sum of one-third of the masses of tetrahedra connected to the gridpoint. If the body is at equilibrium, or in steady-state flow (e.g., plastic flow), \(\sum\ F_i\) on the node will be zero. Otherwise, the node will be accelerated according to the finite difference form of Newton’s second law of motion:
where the superscripts denote the time at which the corresponding variable is evaluated.
During each timestep, strains and rotations are related to nodal displacements in the usual fashion:
Notice that, due to the incremental treatment, Equation (58) does not imply a restriction to small strains.
The constitutive relations for deformable blocks are used in an incremental form so that implementation on nonlinear problems can be accomplished easily. The actual form of the equations is
where:
\(\lambda\), \(\mu\) are the Lame constants;
\(\Delta \sigma_{ij}^e\) are the elastic increments of the stress tensor;
\(\Delta \epsilon_{ij}\) are the incremental strains;
\(\Delta \epsilon_v = \Delta \epsilon_{11} + \Delta \epsilon_{22}\) is the increment of volumetric strain; and
\(\delta_{ij}\) is the Kronecker delta function.
Nonlinear and post-peak strength models are readily incorporated into the code in a direct way without recourse to devices such as equivalent stiffnesses or initial strains, which need to be introduced into matrix-oriented programs to preserve linearity dictated by the matrix formulation. In an explicit program, however, the process is much simpler: After each timestep, the strain state of each zone is known. The program then needs to know the stress in each zone in order to proceed to the next timestep. The stress is uniquely defined by the stress-strain model, whether it is a linearly elastic relation or a complex, nonlinear and post-peak strength model.
The basic failure model for blocks in 3DEC is the Mohr-Coulomb failure criterion with a nonassociated flow rule. Other nonlinear plasticity models are also available (see the i Constitutive Models section for a description of the block material models).
Accurate Modeling of Plastic Collapse — 3DEC is primarily intended to simulate mechanisms related to movement along discrete features (such as joints and faults) within a rock mass. However, in many problems, the failure and collapse of intact material (for example, roof collapse or sloughing of sidewalls of excavations) must also be accommodated in the model.
When using the block plasticity models, described in the i Constitutive Models section, it is important to recognize that an overestimation of the collapse load may be calculated for the constant-strain tetrahedral elements in 3DEC. A common problem that occurs in modeling of materials undergoing active collapse is the incompressibility condition of plastic flow. This condition is sometimes referred to as “mesh-locking” or “excessively stiff” elements, and is discussed in detail by Nagtegaal et al. (1974). The problem arises as a condition of local mesh incompressibility, which must be satisfied during flow, resulting in over-constrained elements.
One method to overcome this problem is referred to as “mixed discretization” (see Marti and Cundall 1982). The principle of the mixed discretization technique is to give the element more volumetric flexibility by proper adjustment of the first invariant of the tetrahedra strain-rate tensor.[2] (This invariant gives a measure of the rate of dilation of the constant strain-rate tetrahedron.) In the approach, a coarser discretization in zones is superposed to the tetrahedral discretization, and the first strain-rate invariant of a particular tetrahedron in a zone is evaluated as the volumetric-average value over all tetrahedra in the zone. The method is illustrated in Figure 15. In the particular mode of deformation sketched there, individual constant strain-rate elements will experience a volume change incompatible with a theory of incompressible plastic flow. In this example, however, the volume of the assembly of tetrahedra (i.e., the zone) remains constant, and application of the mixed discretization process allows each individual tetrahedron to reflect this property of the zone, hence reconciling its behavior with that predicted by the theory.
Figure 15: Mode of deformation for which mixed discretization would be most efficient.
Mixed discretization is available in 3DEC, but only for six-sided polyhedra. The procedure cannot be readily adapted for discretization of arbitrarily shaped blocks. (See the block zone generate hexahedra command.)
In 3DEC, a mixed-discretization (m-d) zone corresponds to an assembly of \(n_t\) tetrahedra, as illustrated in Figure 16 for the case \(n_t\) = 5. Consider a particular m-d zone: the strain-rate tensor of a tetrahedron locally labeled \(l\) in that zone is first estimated and then decomposed into deviatoric and volumetric parts:
where \(\eta^{[l]}\) is the deviatoric strain-rate tensor, and \(\xi^{[l]}\) is the strain-rate first invariant,
Figure 16: An 8-node zone with 2 overlays of 5 tetrahedra in each overlay.
The first invariant for the zone is then calculated as the volumetric average value of the first invariant over all tetrahedra in the zone:
where \(V^{[k]}\) is the volume of the tetrahedron, \(k\). Finally, the tetrahedron strain-rate tensor components are calculated from
Dilatant constitutive laws will produce changes in mean normal stress when yielding occurs. For a consistent technique, the first invariant of the stress tensor, derived after application of the strain-rate increment, must also be evaluated as a volumetric average for the zone. In this process, the stress tensor of a particular tetrahedron, \(l\), in a zone is first estimated and decomposed into deviatoric and volumetric parts,
where \([s]^{[l]}\) is the deviatoric strain-rate tensor, and \({\sigma}^{[l]}\) is the mean normal stress,
The first invariant for the zone is calculated as the volumetric average value over all tetrahedra in the zone:
Finally, the tetrahedron stress-rate tensor components are calculated using
The calculation of nodal forces (based on evaluation of strain rates and stresses) is carried out using a combination of the two overlays. The advantage of the two-overlay approach is to ensure symmetric zone response for symmetric loading. Mixed discretization is carried out over the combination of two overlays, and nodal force computations are evaluated by averaging over the two overlays.
Nodal Mixed Discretization for a Tetrahedral Grid
One of the difficulties associated with low-order elements (such as constant strain) is that they exhibit “volumetric locking” in the analysis of problems in which a constitutive constraint is imposed on the volumetric behavior of the material. Typical problems in this category involve modeling of plastic flow, with zero dilation. Nagtegaal et al. (1974) have explained that the difficulty arises because, under the incompressibility condition, certain classes of meshes are over-constrained. In the example of the Prandtl wedge problem, the use of constant strain elements may not only overpredict the bearing load, but may prove to be unreliable in predicting a limit load capacity at all.
One way to resolve the issue is to increase the order of the element. However, a drawback of introducing additional degrees of freedom (in a higher-order element formulation, such as proposed byKey et al. (1999)) is the possible occurrence of hourglass modes of deformation (combinations of nodal displacements that produce no strain, and are thus unopposed by stresses). These modes are not observed in reality, and application of complicated correction terms in the element formulation is often needed to prevent them from occurring. Another drawback of the higher-order element approach is the increased complexity of the algorithm for contact detection, boundary recognition for application of boundary conditions, and model formulation in general.
To circumvent these difficulties, Marti and Cundall (1982) have proposed a procedure that reduces the number of constraints on plastic flow, and at the same time keeps the order of the element low, thus preventing unwanted hourglassing. The technique is called “mixed discretization” because the discretization for the isotropic and deviatoric parts of the stress and strain tensors are different. In essence, the deviatoric behavior is defined on an element basis (triangle or tetrahedra), while the volumetric behavior is averaged over an assembly of elements, referred to as zones (quadrilatera or hexahedra).
Description of the technique
Nodal mixed discretization (NMD) is a variation on the mixed discretization scheme, in which averaging of the volumetric behavior is carried out on a node, rather than a zone, basis. The procedure is applied on the triangle- or tetrahedral-based mesh; it does not require the assembly of elements into zones. Also, the constitutive model is called on an element basis, as usual (no call is made on a node basis for the volumetric behavior, as in the average nodal pressure (ANP) formulation described by Bonet and Burton (1998)). The procedure involves a nodal mixed discretization on strain, and one on stress. Both steps are considered below.
First we recall the general calculation sequence embodied in FLAC and 3DEC:
- Nodal forces are calculated from stresses, applied loads and body forces (velocity and displacement vary linearly; stress and strain are constant within an element).
- The equations of motion are invoked to derive new nodal velocities and displacements.
- Element strain rates are derived from nodal velocities.
- New stresses are derived from strain rates, using the material constitutive law.
In the NMD technique, the calculation sequence is respected. However, an averaging procedure is carried out on strain rates (end of step 3) and on stress increments (end of step 4), as described below.
Nodal mixed discretization on strain
The strain rate \(\dot \varepsilon_{ij}\) is derived from nodal velocities, as usual. The strain rate is then partitioned into deviatoric, \(\dot e_{ij}\), and volumetric, \(\dot e\), components:
where \(\delta_{ij}\) is the Kroenecker delta.
A nodal volumetric strain rate (defined as the weighted average of the surrounding element values) is calculated using the formula
where \(m_n\) are the elements surrounding node \(n\), and \(V_e\) is the volume of element \(e\).
After nodal volumetric strain rate values are obtained, a mean value for the element \(\bar{\dot e}\) is calculated by taking the average of nodal values:
where \(d\) = 3 for a triangle and 4 for a tetrahedral.
Finally, the element strain rate is redefined by superposition of the deviatoric part and volumetric average:
The constitutive model is called to derive new stresses (from strain rates) and previous stresses.
Nodal mixed discretization on stress
Consider an incremental volumetric constitutive stress-strain law, which, for small increments, can be linearized in the form
where \(\dot e^p\) stands for plastic, volumetric strain increment, and the value is nonzero for dilatant/contractant material. The associated nodal forces must be consistent with the assumptions taken to define the element kinematics. To enforce this, a nodal mixed discretization procedure is applied on the term \(K\dot e^p\), as described below. For convenience, we refer to the term \(K\dot e^p\) as \(\dot \sigma^p\); With this convention, Equation (72) may be expressed as
The value \(\dot \sigma^p\) is a standard quantity evaluated in the constitutive model procedure.
The technique for nodal mixed discretization on stress is similar to the one applied for strain. First, nodal values for \(\dot \sigma^p\) are calculated as the weighted average of the surrounding element values,
After the nodal values \(\dot \sigma_n^p\) are obtained, a mean value for the element \(\bar{\dot \sigma^p}\) is calculated by taking the average of nodal values,
where, again, \(d\) = 3 for a triangle and 4 for a tetrahedra.
Finally, the stresses calculated by the constitutive model are corrected by substituting \(\bar{\dot \sigma^p}\) for \(\dot \sigma^p\) in the returned value:
Clearly, the nodal mixed discretization on stress will only be relevant for dilatant/contractant materials.
It is important to note that the averaging manipulations involved in the NMD technique are done independent of the constitutive model formulation, which remains unaffected. (The plastic volumetric “stress correction” \(\dot \sigma^p\), calculated by the constitutive model, is simply passed back from the constitutive model as one of the state variables.)
Also, it may be shown that, because of the way the averaging procedure is carried out, the resulting tetrahedral formulation gives the correct behavior in patch test situations where \(\dot \sigma\) is uniform, and thus equal, for all elements.
The nodal averaging procedure involved in the NMD technique removes the excessively constrained kinematics of the linear velocity element.
Higher-Order Tetrahedral Elements
The standard 3DEC zones are 4-node tetrahedra, assuming linear displacement interpolation functions. The command block zone generate high-order-tetra transforms a mesh of standard zones into a mesh of higher order tetrahedra, with 10 nodes, based on quadratic displacement interpolations functions. For this purpose, new nodes are created at the midpoint of every zone edge. The higher-order element formulation allows a quadratic displacement field to be represented inside the zone and also on the zone faces. However, for purposes of contact calculations, and for plotting, the block boundary is approximated by a mesh of triangular faces. Each face of a 10-node tetrahedron is divided into 4 plane triangles.
The higher-order tetrahedral zones are based on a standard finite element formulation, similar to the one used for the 20-node brick elements. For the 10-node tetrahedra, there are 4 Gauss points, where strains and stresses are calculated. The element is formulated to be compatible with the 3DEC explicit solution method. Therefore, nodal forces are obtained from the Gauss point stresses by numerical integration, so that the stiffness matrix is not calculated. At each Gauss point, the constitutive model adopted for the block is applied at every step to obtain the new stresses, as is done for the standard zones.
For printing and plotting purposes, a mesh of 4 standard zones is superimposed on each 10-node tetrahedron. The Gauss point stresses are transferred to these zones, so that the standard printing contouring commands in 3DEC may still be used.
The higher-order zones provide a better stress approximation with coarser meshes. They also provide a more accurate solution for block plasticity problems. These elements are therefore suited for problems in which the stress distribution inside the blocks is important, or significant block yielding takes place. For typical blocky systems, where the contact behavior is dominant, the higher-order zones may be less appropriate, given the approximations involved in replacing the curved, deformed block faces by plane triangles, as assumed in the contact calculations.
Mechanical Damping
Mechanical damping is used in the distinct element method to solve two general classes of problems: static (non-inertial) solutions and dynamic solutions. A different form of damping is used for each class. For static analysis, the approach is conceptually similar to dynamic relaxation, proposed by Otter et al. (1966). The equations of motion are damped to reach a force equilibrium state as quickly as possible under the applied initial and boundary conditions. Damping is velocity-proportional (i.e., the magnitude of the damping force is proportional to the velocity of the blocks).
The use of velocity-proportional damping in standard dynamic relaxation involves three main difficulties:
- The damping introduces body forces, which are erroneous in “flowing” regions and may influence the mode of failure in some cases.
- The optimum proportionality constant depends on the eigenvalues of the matrix, which are unknown unless a complete modal analysis is done. In a linear problem, this analysis needs almost as much computer effort as the dynamic relaxation calculation itself. In a nonlinear problem, eigenvalues may be undefined.
- In its standard form, velocity-proportional damping is applied equally to all nodes (i.e., a single damping constant is chosen for the whole model). In many cases, a variety of behavior may be observed in different parts of the model; for example, one region may be failing while another is stable. For these problems, different amounts of damping are appropriate for different regions.
In an effort to overcome one or more of these difficulties, alternative forms of damping may be proposed. In soil and rock, natural damping is hysteretic; if the slope of the unloading curve is higher than that of the loading curve, energy may be lost. This type of damping can be produced numerically, but there are at least two difficulties. First, the precise nature of the hysteretic curve is often unknown for complex loading-unloading paths. This is particularly true for soils, which are typically tested with sinusoidal stress histories. Cundall (1976) reports that very different results are obtained when the same energy loss is accounted for by different types of hysteretic loops. Second, “ratcheting” can occur (i.e., each cycle in the oscillation of a body causes irreversible strain to be accumulated). This type of damping has been avoided, since it increases path-dependence and makes the results more difficult to interpret.
Two alternative forms of velocity-proportional damping are provided in 3DEC. The first is a numerical servo-mechanism, termed adaptive global damping, and is described by Cundall (1982). Adaptive global damping is used to adjust the damping constant automatically. Viscous damping forces are used, but the viscosity constant is continuously adjusted in such a way that the power absorbed by damping is a constant proportion of the rate of change of kinetic energy in the system. The adjustment to the viscosity constant is made by a numerical servo-mechanism that seeks to keep the following ratio, \(R\), equal to a given ratio (e.g., 0.5):
where:
\(P\) is the damping power for a node;
\(\dot E_k\) is the rate of change of nodal kinetic energy; and
\(\sum\) represents the summation over all nodes.
This form of damping overcomes difficulty (2) above, and partially overcomes (1), since, as a system approaches steady state (equilibrium or steady flow), the rate of change of kinetic energy approaches zero, and consequently, the damping power tends toward zero (see Cundall 1982).
In order to overcome all three difficulties, 3DEC provides another form of damping, in which the damping force on a node is proportional to the magnitude of the unbalanced force. For this scheme, referred to as local damping, the direction of the damping force is such that energy is always dissipated. For deformable blocks, the equation of motion Equation (57) is replaced by an equation that incorporates local damping:
where
where \(\alpha\) is a constant (set to 0.8 in 3DEC ), and \(m_n\) is the nodal mass. A similar equation is used in place of Equation (49) to apply local damping to translational and angular velocities of rigid blocks.
This type of damping is equivalent to a local form of adaptive damping. In principle, the three difficulties reported above are addressed: body forces vanish for steady-state conditions; the magnitude of damping constant is dimensionless and is independent of properties or boundary conditions; and the amount of damping varies from point to point (Cundall 1987, pp. 134-135).
UDEC calculations using local damping and adaptive global damping are compared by Cundall (1987); the methods are shown to converge to the same solution. Local damping may be preferred for analyses involving sudden load changes or progressive failure (such as caving of many blocks), for which a different amount of damping is required in different regions of the model. Analyses with local damping are observed to be slightly underdamped in general. Adaptive global damping is the default damping mode for static analysis in 3DEC. Adaptive global damping is assigned with the damping auto command, and local damping is assigned with the damping local command.
Combined Damping — A variation on local damping is also provided, for situations in which the steady-state solution includes a significant uniform motion, as may occur in a creep simulation. This damping is called combined damping. For this special case, combined damping removes kinetic energy more efficiently than local damping.
The damping formulation described by Equation (78) is only activated when the velocity component changes sign. In situations where there is significant uniform motion (in comparison to the magnitude of oscillations that are to be damped), there may be no “zero-crossings,” and hence no energy dissipation.
In order to develop a damping formulation that is insensitive to rigid-body motion, consider periodic motion superimposed on steady motion:
where \(V\) is the maximum periodic velocity, \(\omega\) is the angular frequency and \(\dot u_\circ\) is the superimposed steady velocity. Differentiating twice, and noting that \(m\ddot u = F\),
In Equation (81), \(\dot F\) is proportional to the periodic part of \(\dot u\), without the constant \(\dot u_\circ\). We may substitute \(-{\rm sgn}(\dot F)\) for the damping force in Equation (78) to obtain the same damping force, if the motion is periodic:
This equation is insensitive to a constant offset in velocity, since \(\dot F\) does not involve \(\dot u_\circ\). In practice, Equation (82) is not as efficient as the local damping force term, (79), if the motion is not strictly periodic. However, the combination of both formulas in equal proportions gives good results:
This form of damping should be used if there is significant rigid-body motion in a system, in addition to oscillatory motion to be dissipated. For this reason, combined damping is the default damping mode for creep analysis. In most cases, local damping is preferred because combined damping is found to dissipate energy at a slower rate than local damping, based on velocity.
For a dynamic analysis, the damping in the numerical simulation should approximately reproduce the energy losses in the natural system when subjected to a dynamic loading. As mentioned above, in soil and rock, natural damping is mainly hysteretic (i.e., independent of frequency). It is difficult to reproduce this type of damping numerically because of the problem with path-dependence, as described previously. As an alternative, Rayleigh damping is used in 3DEC. This method of damping for dynamic analysis is described in the Dynamic Analysis section.
Numerical Stability
The solution scheme used for the distinct element method is conditionally stable. A limiting timestep that satisfies the stability criterion for both the calculation of internal block deformation and inter-block relative displacement is determined. The timestep required for the stability of block deformation computations is estimated as
where:
\(m_i\) is the mass associated with block node \(i\); and
\(k_i\) is the measure of stiffness of the elements surrounding the node.
The ratio of mass to stiffness is related to the highest eigenfrequency, \(\omega_{\rm max}\), of a linear elastic system.
The stiffness term, \(k_i\), must account for both the stiffness of the intact rock and that of the discontinuities. It is calculated as the sum of the two components:
The first term on the right-hand side represents the sum of the contributions of the stiffness of all elements connected to node \(i\), which are estimated as
where:
\(K\) and \(G\) are the bulk and shear elastic moduli of the block material, respectively;
\(b_{\max}\) is the largest zone edge; and
\(h_{\min}\) is the minimum height of the tetrahedral element.
The joint stiffness term, \(k_ji\), exists only for nodes located on the block boundary, and is taken as the product of the normal or shear joint stiffnesses (whichever is larger) and the sum of the areas of the two block face segments adjacent to node \(i\).
For calculations of inter-block relative displacement, the limiting timestep is calculated, by analogy to a simple degree-of-freedom system, as
where:
\(M_{\min}\) is the mass of the smallest block in the system; and
\(K_{\max}\) is the maximum contact stiffness.
The term frac is a user-supplied value that accounts for the fact that a single block may be in contact with several blocks simultaneously. A typical value for frac is 0.1.
The controlling timestep for a distinct element analysis is
Mass (Density) Scaling
Some way of increasing the timestep is desirable, in order to reduce computer time, even though explicit calculations execute very rapidly per timestep. One way to do this is by scaling the mass (or density) of the solid material. It may be noted that the value of inertial density is irrelevant to the modeling of static systems, provided that gravity forces are correctly preserved. For nearly static systems (i.e., ones which only evolve slowly with time), inertial densities may be increased until they begin to become appreciable compared to other forces in the system. The system response will not be significantly modified if inertial forces are low. The reason for wanting to increase the density is that the critical timestep, \(\Delta t\), may also be increased because it is determined by density, \(\rho\):
This procedure is called density scaling. Density scaling is only effective at improving convergence if the model is nonuniform (i.e., if the natural timesteps differ for different parts of the model). Changing the density in a uniform model does nothing to improve convergence.
Density scaling is the default mode for static analysis with 3DEC. It may also be selected and adjusted by the user with the block mechanical mass-scale command. Mass scaling is turned on automatically when local damping is specified (via block mechanical damping local) or when adaptive global damping is specified (via block mechanical damping global). For most problems, a scale factor based on the average block mass or zone mass in the model provides the most rapid convergence.
Partial Mass Density Scaling for Dynamic Analysis
Density scaling is a technique (used in 3DEC in quasi-static calculations) that substantially improves the efficiency in obtaining solutions to large problems. As in quasi-static problems, inertial forces are not important; the gridpoint masses can be scaled for optimal numerical convergence without affecting the solution. In dynamic analyses, however, global scaling cannot be used. For complex jointed systems, very small block zones, which require very small timesteps for numerical stability of the explicit algorithm, are created during the automatic meshing procedure. This makes some dynamic solutions extremely time-consuming. However, as these zones may be very small, with very small masses, it is possible to introduce some density scaling only for these zones in such a way that the change of the system inertia is negligible. This scheme of partial density scaling was implemented in 3DEC in such a way that the user controls the amount of scaling to be introduced. Given the timestep calculated by the code, the user specifies the desired timestep with the command block mechanical mass-scale timestep f. JH: is this the correct substitution (old docs had “mscale part dt”)? And see single keyword reference in next paragraph – fix as needed then delete this comment please
This command specifies that only the amount of density scaling required to achieve the timestep specified by timestep f is to be applied to the system. When a model cycle command is given, a message indicating the number of gridpoint masses that were scaled and the amount of additional mass introduced is printed.
References
Bardet, J.-P., and R. F. Scott. “Seismic Stability of Fractured Rock Masses with the Distinct Element Method,” in Research and Engineering Applications in Rock Masses (Proceedings of the 26th U.S. Symposium on Rock Mechanics), Vol. 2, pp. 139-150. Boston: A. A. Balkema (1985).
Bonet, J., and A. J. Burton. “A Simple Averaged Nodal Pressure Tetrahedral Element for Nearly Incompressible Dynamic Explicit Applications,” Commun. Numer. Meth. Engng., 14, 437-449 (1998).
Butkovich, T. R., O. R. Walton and F. E. Heuze. “Insights in Cratering Phenomenology Provided by Discrete Element Modeling,” in Key Questions in Rock Mechanics: Proceedings of the 29th U.S. Symposium (University of Minnesota, June 1988), pp. 359-368. Rotterdam: A. A. Balkema (1988).
Cundall, P. A. “A Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock Systems,” in Proceedings of the Symposium of the International Society of Rock Mechanics (Nancy, France, 1971), Vol. 1, Paper No. II-8 (1971).
Cundall, P. A. “Adaptive Density-Scaling for Time-Explicit Calculations,” in Proceedings of the 4th International Conference on Numerical Methods in Geomechanics (Edmonton, Canada, 1982), pp. 23-26 (1982).
Cundall, P. A. “Distinct Element Models of Rock and Soil Structure,” in Analytical and Computational Methods in Engineering Rock Mechanics, Chapter 4, pp. 129-163. E. T. Brown, ed. London: George Allen and Unwin. (1987).
Cundall, P. A. “Explicit Finite Difference Methods in Geomechanics,” in Numerical Methods in Engineering (Proceedings of the EF Conference on Numerical Methods in Geomechanics (Blacksburg, Virginia, June 1976), Vol. 1, pp. 132-150 (1976).
Cundall, P. A. “Formulation of a Three-Dimensional Distinct Element Model – Part I: A Scheme to Detect and Represent Contacts in a System Composed of Many Polyhedral Blocks,” Int. J. Rock Mech., Min. Sci. & Geomech. Abstr., 25, 107-116 (1988).
Cundall, P. A. “Rational Design of Tunnel Supports: A Computer Model for Rock Mass Behaviour Using Interactive Graphics for the Input and Output of Geometrical Data,” U. S. Army Corps of Engineers, Missouri River Division, Technical Report MRD-2-74; NTIS Report No. AD/A-001 602 (1974).
Cundall P. A. “UDEC – A Generalized Distinct Element Program for Modelling Jointed Rock,” Report PCAR-1-80, Peter Cundall Associates Report, European Research Office, U.S. Army, Contract DAJA37-79-C-0548 (1980).
Cundall, P. A., and R. D. Hart. “Development of Generalized 2-D and 3-D Distinct Element Programs for Modeling Jointed Rock,” Itasca Consulting Group; Misc. Paper SL-85-1, U.S. Army Corps of Engineers (1985).
Cundall, P. A., and R. D. Hart. “Numerical Modeling of Discontinua,” Engr. Comp., 9(2), 101-113 (1992).
Cundall, P. A., et al. “Computer Modeling of Jointed Rock Masses,” U.S. Army EngineerWaterways Experiment Station, Vicksburg, Mississippi, Tech. Report N-78-4 (August 1978).
Cundall P. A., and O. D. L. Strack. “A Discrete Numerical Model for Granular Assemblies,” Geotechnique, 29, 47-65 (1979b).
Cundall, P. A., and O. D. L. Strack. “Modeling of Microscopic Mechanisms in Granular Material,” in Mechanics of Granular Materials: New Models and Constitutive Relations, pp. 137-149. Amsterdam: Elsevier Scientific Publications, B.V. (1983).
Cundall P. A., and O. D. L. Strack. “The Distinct Element Method as a Tool for Research in Granular Media, Part I,” University of Minnesota, Department of Civil and Mineral Engineering, Report to National Science Foundation, Grant ENG 76-20711 (1979a).
Detournay, C., and E. Dzik. “Nodal Mixed Discretization for Tetrahedral Elements” in FLAC and Numerical Modeling in Geomechanics (Proceedings of the 4th International FLAC Symposium, Madrid, Spain, May 29-31 2006). P. Varona & R. Hart, eds. Minneapolis: Itasca, 343-350 (2006).
Goodman R. E., and G. Shi. Block Theory and Its Application to Rock Engineering. New Jersey: Prentice-Hall (1985).
Hahn, J. K. “Realistic Animation of Rigid Bodies,” Computer Graphics, 24(4), 299-308 (1988).
Hart, R., P. Cundall and J. Lemos. “Formulation of a Three-Dimensional Distinct Element Model – Part II: Mechanical Calculations for Motion and Interaction of a System Composed of Many Polyhedral Blocks,” Int. J. Rock Mech., Min. Sci. & Geomech. Abstr., 25, 117-126 (1988).
Hart, R. D., J. Lemos and P. Cundall. “Block Motion Research: Analysis with the Distinct Element Method,” Itasca Consulting Group/Agbabian Associates, DNA-TR-88-34-V2 (December 1987).
Hart, R. D., et al. “Tunnel Prediction Using Distinct Elements: Volume II – Computer Code Modification and Verification,” DNA, Technical Report DNA-TR-90-56-V2 (December 1990).
Heuzé, F. E., et al. “Analysis of Explosions in Hard Rock: The Power of Discrete Element Modeling,” Lawrence Livermore Laboratory, Preprint UCRL-JC-103498 (March 1990).
Hocking, G., G. G. W. Mustoe and J. R. Williams. CICE Discrete Element Code – Theoretical Manual. Lakewood, Colorado: Applied Mechanics Inc. (1985).
Itasca Consulting Group Inc. Particle Flow Code in 2 Dimensions and Particle Flow Code in 3 Dimensions, Version 4.0. Minneapolis: ICG (2008).
Itasca Consulting Group Inc. Universal Distinct Element Code, Version 5.0. Minneapolis: ICG (2011).
Key, S. W. “A Data Structure for Three-Dimensional Sliding Interfaces,” in Computational Mechanics ’86: Theory and Applications (Proceedings of the International Conference on Computational Mechanics, Tokyo, Japan, 1986), pp. V1-211 to V1-216. Tokyo: Springer-Verlag (1986).
Key, S. W., et al. “A Suitable Low-Order, Tetrahedral Finite Element for Solids,” Int. J. Numer. Meth. Engng., 44, 1785-1805 (1999).
Lemos, J. “A Distinct Element Model for Dynamic Analysis of Jointed Rock with Application to Dam Foundations and Fault Motion.” Ph.D. Thesis, University of Minnesota (1987).
Lemos, J. V., R. D. Hart and P. A. Cundall. “A Generalized Distinct Element Program for Modeling Jointed Rock Mass (A Keynote Lecture),” in Proceedings of the International Symposium on Fundamentals of Rock Joints (Björkliden, Sweden, September 1985), pp. 335-343. Luleå, Sweden: Centek Publishers (1985).
Lorig, L. J., and P. A. Cundall. “Modeling of Reinforced Concrete Using the Distinct Element Method,” in Fracture of Concrete and Rock, pp. 459-471. Bethel, Conn.: SEM (1987).
Marti, J., and P. A. Cundall. “Mixed Discretization Procedure for Accurate Solution of Plasticity Problems,” Int. J. Num. Methods & Analy. Methods in Geomech., 6, 129-139 (1982).
Nagtegaal, J. C., D. M. Parks and J. R. Rice. “On Numerically Accurate Finite Element Solutions in the Fully Plastic Range,” Comp. Meth. Appl. Mech., 4, 153-177 (1974).
Otter, J. R. H., A. C. Cassell and R. E. Hobbs. “Dynamic Relaxation (Paper No. 6986),” Proc. Instn. Civ. Engrs., 35, 633-656 (1966).
Papadopoulos, J. M. “Incremental Deformation of an Irregular Assembly of Particles in Compressive Contact.” Ph.D. Thesis, M.I.T, Department of Mechanical Engineering (1986).
Plesha, M. E., and E. C. Aifantis. “On the Modeling of Rocks with Microstructure,” in Rock Mechanics – Theory-Experiment-Practice (Proceedings of the 24th U.S. Symposium on Rock Mechanics, Texas A&M University, 1983), pp. 27-35. New York: Association of Engineering Geologists (1983).
Potyondy, D. O., P. A. Cundall and C. Lee. “Modeling Rock Using Bonded Assemblies of Circular Particles,” in Rock Mechanics Tools and Techniques, pp. 1937-1944. Rotterdam: A. A. Balkema (1996).
Shi, G.-H. “Discontinuous Deformation Analysis – A New Numerical Model for the Statics and Dynamics of Block Systems,” Lawrence Berkeley Laboratory, Report to DOE OWTD, Contract AC03-76SF0098 (September 1988); also Ph.D. Thesis, University of California, Berkeley (August 1989).
Walton, O. R. “Particle Dynamic Modeling of Geological Materials,” Lawrence Livermore National Laboratory, Report UCRL-52915 (1980).
Walton, O. R., et al. “Particle-Dynamics Calculations of Gravity Flows of Inelastic, Frictional Spheres,” in Micromechanics of Granular Material, pp. 153-161. Amsterdam: Elsevier Science Publishers (1988).
Warburton, P. M. “Vector Stability Analysis of an Arbitrary Polyhedral Rock Block with any Number of Free Faces,” Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 18, 415-427 (1981).
Williams, J. R., and G. G. W. Mustoe. “Modal Methods for the Analysis of Discrete Systems,” Computers & Geotechnics, 4, 1-19 (1987).
Endnotes
| [1] | The term “hourglassing” comes from the shape of the deformation pattern of elements within a mesh. For polyhedra with more than four nodes, combinations of nodal displacements which produce no strain and result in no opposing forces exist. The resulting effect is unopposed deformations of alternating direction. |
| [2] | This invariant gives a measure of the rate of dilation of the constant strain-rate tetrahedron. |
| Was this helpful? ... | PFC © 2021, Itasca | Updated: Feb 25, 2024 |