Research Article  Open Access
Discontinuous Deformation Analysis Enriched by the Bonding Block Model
Abstract
The discontinuous deformation analysis (DDA) has been extensively applied in geotechnical engineering owing to its salient merits in the modeling of discontinuities. However, this method assumes a constant stress field within every block and hence cannot provide reliable estimation for block deformations and stresses. This paper proposes a novel scheme to improve the accuracy of the DDA. In our method, advanced subdivision is introduced to represent a block as an assembly of triangular or quadrilateral elements, in which overlapped element edges are separated from each other and are glued together by bonding springs. The accuracy and the effectiveness of the proposed method are illustrated by three numerical experiments for both continuous and discontinuous problems.
1. Introduction
The discontinuous deformation analysis (DDA) [1, 2] is developed for the modeling of the statics and dynamics problems in geological engineering. This method represents a jointed rock mass as an assembly of blocks with constantly changing deformations and contact status. At every contact point, a normal spring and a shear spring are applied. Such contact springs must satisfy no tension in an open contact and no penetration in a close contact, which are fulfilled through repeated openclose iterations. In the literature [3–8], the augmented Lagrangian method is also introduced into the DDA to increase the accuracy for contact computation. The DDA uses an incremental procedure to solve block movements and deformations on a given time mesh. In every time step, the simultaneous equilibrium equations are formed by the submatrices derived from minimizing all sources of potential energies. Due to the above unique merits, this method has been extensively applied in the analysis of seismic landslides [9–11], crack propagations [12–14], and hydraulic fractures [15, 16].
However, the DDA employs the firstorder polynomial to approximate block displacements, and hence the stress within every block keeps invariable. Such assumption leads to the difficulty to implement this method to simulate block cracks, which require accurate stress estimation. In order to improve the accuracy of the DDA, a multitude of methods have been developed. A direct way is to use second or higherorder polynomial functions to approximate block displacements [17–22]. But the accuracy of this method is only adequate for regularly shaped blocks. Another way is the socalled subblock DDA [3, 13, 16, 23], which subdivides blocks into smaller ones by preassumed artificial joints and applies contact springs along the interfaces to prevent the relative movements. Since contact detections and openclose iterations are also required for the contacts between subblocks, the computational burden of this method is very heavy, especially when adopting finer subdivisions.
The accuracy of the DDA can also be improved by coupling with the FEM. In the literature [24, 25], the DDA and the FEM codes are integrated into a computer program to alternately figure out block deformations and contact forces. Besides, a FEM postprocessing technique is introduced to recover block stresses from DDA solution in [14, 26]. A disadvantage of these two methods lies in that, in every time step, both DDA equations and FEM/postprocessing equations have to be solved, and their results also need to be passed to each other. The additional calculation tasks reduce the computing efficiency. Different from the above coupling schemes, the method presented in [27] divides a problem zone into two domains, on which the DDA and the FEM are implemented, respectively. Along the domain interface, the socalled line blocks are introduced to transfer the contact forces computed by the DDA and the deformation of FEM domain. However, this method can only analyze statics problems and need complicated displacement control algorithm to ensure the convergence of the computation. Alternatively, the nodalbased DDA [28–33] incorporates finite element type mesh on every block and uses nodal displacements as the unknowns to be determined. But there still exists the difficulty due to the constraint on nodal displacements which must continue across element interfaces, especially in large deformation problems.
In this paper, we propose a novel scheme referred to hereafter as the bonding block model (BBM) in the framework of the DDA to improve the accuracy of block deformations and stresses. Compared with the standard DDA, our method represents every nonrigid block as an assembly of triangular or quadrilateral elements, as shown in Figure 1, where overlapped element edges are separated from each other and are glued together by the introduced bonding springs. Contacts are allowed to occur along the element edges that compose block boundaries. At every contact point, a normal contact spring and a shear contact spring are applied, and our method retains the openclose iteration procedure developed in the standard DDA to handle their installation. Benchmark numerical experiments involving both continuous and discontinuous problems are conducted. The results are compared with analytical solutions and ANSYS simulations to illustrate the accuracy and effectiveness of the proposed method.
Although both of our method and the subblock DDA subdivide blocks into smaller partitions, the main difference lies in that we connect the elements within a block by bonding springs rather than the contact springs applied on subblock interfaces. In this way, it is avoided to implement contact detections and time consuming openclose iterations along artificial joints. Compared with the coupling DDAFEM schemes, our method is a dynamics one, and, in every time step, block deformations and contact effects can be simultaneously figured out without alternately executing the DDA subroutine and the FEM/postprocessing subroutine. The proposed method is also distinguished from the nodalbased DDA by the unique discretization, in which overlapped element edges are separated from each other but are not sharing a common grid line, and adjacent elements are connected by bonding springs instead of continuous element shape functions.
The rest of the paper is organized as follows. In the next section, basic concepts of the proposed method are introduced. In Section 3, detailed formulations are derived. In Section 4, three numerical experiments are conducted to illustrate the accuracy and the efficiency of our method for the analyses of both continuous and discontinuous problems. Section 5 concludes this paper.
2. Concepts of the BBM
2.1. Subdivision
For a jointed rock mass, assume it contains blocks occupying the volumes . As shown in Figure 1, the BBM implements an advanced subdivision over blocks to enhance their flexibility. The subdivision process can be schematized from Figure 2(a) to Figure 2(b), where the block is partitioned into the triangular elements , . As illustrated in Figure 2(c), the overlapped element edges are separated from each other. Proceed such subdivision over all blocks, and assume elements are generated, including . On each element, say , the displacement at any point can be approximated by the following function:where is the displacement mode matrix:in which denotes the barycenter coordinate of the element ; is composed of the unknowns to be determined of the element; that is,where , are the translational displacements of the barycenter, is the rotation of the centroid, and , , and denote the strain components of .
(a)
(b)
(c)
Since the approximation introduced above is independent of mesh nodes, arbitrarily shaped elements can be adopted, such as triangles, quadrilaterals, and even original blocks without subdivision. Although the formulations in Section 3 are derived with the illustration figures for triangular elements, they can be directly employed in the analysis by using other element shape options.
2.2. Element Connections
For an element system obtained through the process introduced above, the element edges can be classified into two types: the “boundary edges” that are parts of block boundaries where contacts may occur and the “internal edges” including the element edges inside original blocks. As shown in Figure 2, and are two adjacent elements of the block . From Figure 3(a), it can be recognized that and are boundary edges, while , , , and are internal edges. Just like the DDA, the BBM applies a normal contact spring and a shear contact spring at every contact point and controls their installation through openclose iterations.
(a)
(b)
In order to preserve the displacement consistency between each overlapped internal edge, a couple of bonding springs are applied between the coincided vertexes. As illustrated in Figure 3(b), the overlapped edges and are glued together by two bonding springs applied between and and between and , respectively. The bonding springs between every two elements have the same stiffness. And we use the following formulation to determine the bonding spring stiffness between the elements and :where is a constant, is the average shear modulus of the two elements, and and denote the areas of and , respectively. By using to denote the ratio of the maximum element area with the minimum element area for a subdivision, we can take the value of ranging from 10 (for ) to 1000 (for ).
The stiffness given by (4) is the minimum requirement for the proper implementation of bonding springs. Although larger stiffness is also suitable in theory, very large one is detrimental to the computing convergence and even leads to illconditioned global stiffness matrix. For the simplicity, we use instead of to denote bonding spring stiffness in the following sections.
3. Computing Formulations
The BBM implements an incremental iteration process to solve both dynamics and statics problems of the element system on a given time mesh . In particular, time steps used in statics problems represent the loading steps. Below, the th time step indicates the interval . Obviously, it is equivalent to say the value of a variable at the end of the th step and the one at the beginning of the th step.
For the element system generated in the previous section, the total potential energy is the summation of all sources of potential energies, including the strain energies, inertia effects, external loads, volume forces, bonding loadings, and contact forces. Minimizing the potential energy in the th step leads to the global simultaneous equilibrium equations:where the submatrices and are of the sizes and , respectively, while denotes the unknown vector to be determined, . The detailed formulations of these submatrices, and , are derived as follows.
3.1. Element Elastic Submatrices
The strain energy of the th element can be expressed aswhere the vectors , , and denote the strain, stress, and the initial stress in the th time step, respectively. For isotropic linearly elastic materials, and have the following relationship:in which and indicate Young’s modulus and Poisson’s ratio, respectively. The strain can be expressed by the displacement asSubstituting (7) and (8) into (6), we haveMinimizing with respect to leads towhich are added to the stiffness matrix and the right item of (5), respectively.
3.2. Element Inertia Submatrices
The potential energy contributed by the inertia force acting on the element at the moment can be expressed as follows:where and indicate the displacement and the mass per area, respectively, and the multiplication between the parentheses in the integrand of the equation above thereby computes the density of the inertia force on . Substituting (1) into (11) leads toUsing the Newmark time integrating method, we have the following finite difference approximation:where is the th stepsize and indicates the initial velocity of the th time step and can be updated by the iteration formulation below:Substituting (13) into (12), we then haveMinimizing the right items of the equation above with respect to leads towhere the integrals can be solved by the simplex integral method [34]. The submatrices equation (16) and equation (17) are added to the stiffness matrix and the right item of the global simultaneous equilibrium equations in (5), respectively.
3.3. Bonding Submatrices
In order to preserve the displacement consistency along the element interfaces within a block, each pair of initially overlapped element edges is glued together by two bonding springs applied between the coincided vertexes. As shown in Figure 4(a), the points and are a pair of initially overlapped vertexes of the elements and which are within the same block. A bonding spring with the stiffness is applied between and to glue and together. Denote the coordinate of and at the time node by and , respectively.
(a)
(b)
Firstly, we compute the elongation of the bonding spring between and at the time node . As illustrated in Figure 4(b), is actually the length of and can be computed from the following equality:where denotes the vector at , that is,Assume the displacements of and in the th time step are and , respectively. Then, the coordinates of and at can be decomposed into the incremental form below:Substituting (20a) and (20b) into (19) leads toSubstitute (3) and (1) into the equation above, and then we havewhere and . Substituting (22) into (18), then can be expanded as the following summation:
Next, we compute the strain energy of the bonding spring between and in the th time step . can be expressed as follows:Substitute (23) into the expression above, and becomes the following sum:where
Then, minimize with respect to the unknowns. As is known in the th time step, the derivative of is 0. Minimize , , , and with respect to and , and we haveThese submatrices in the equation above are added into the global stiffness matrix in (5). Minimizing and with respect to and leads towhich are added to the right side of (5).
3.4. Bonding Submatrices at Fixed Points
As the boundary conditions, some elements are fixed at specific points which are called fixed points in the original DDA. As shown in Figure 5(a), the point of the element is a fixed point, . Denote the coordinates of at the time node by , . Assume there exists a fictitious element which is immersed in the support wall and always keeps a statics state during a computation. Such fictitious element does not contribute potential energy to the element system. Assume the vertex of overlaps with the fixed point at the initial time . In order to resist the displacement of the point , a bonding spring with the stiffness is applied between and . The stiffness can be determined by the following formulation:where is the same factor used in (4).
(a)
(b)
As shown in Figure 5(b), the bonding spring between and has the elongation at the time node . Denote the coordinate of by , which keeps constant during a computation. Since is actually the length of , we havewhere denotes the vector at ; that is,Assume the displacement of in the th time step is . Then, the coordinates of at can be written as , which is substituted into (31); that is,Substitute (3) and (1) into the equation above, and then we havewhere indicates . Substituting (33) into (30) leads toThen we can express the strain energy of the bonding spring between and aswhereObviously, the minimization of with respect to is zero. Then, we minimize and :where and are added to the stiffness matrix and the right item of (5).
3.5. Line Loading Submatrices
Assume a loading is applied along a line upon the boundary of the element , . The coordinate of any point on the line at the time moment can be determined by the following parametric equations:where and are the coordinates of the ending points of the line, whose length at the end of the last time step can be computed bywhere time step . Then, we can obtain the potential energy contributed by the loading :Substituting (3) and (1) into (40) leads toMinimize with respect to , and we haveThe submatrix is added to the right item of the global equation (5). If , which is constant, (42) has the analytical form:where and is the coordinate of the barycenter of the element .
3.6. Volume Force Submatrix
When the element bears the body loading , , the potential energy due to the force is calculated asSubstituting (1) into the equations above leads toMinimize with respect to , and we haveThe submatrix above is added to the right side of the global simultaneous equation (5). When the body force is constant, say , (46) becomeswhere denotes the area of the element .
3.7. Contact Submatrices
The element edges along the boundaries of different blocks may contact with each other. This kind of interactions is allowed in three modes: vertextovertex, edgetoedge, and vertextoedge contacts. In twodimensional problems, both vertextovertex and edgetoedge contacts can be decomposed into couples of vertextoedge ones. The BBM adopts the “rigid” contact model. In order to prevent penetration and relative sliding along the interfaces of a close contact, a normal contact spring and a shear contact spring have to be applied at the contact point.
During each time step, the simultaneous equilibrium equations in (5) have to be assembled and solved repeatedly whenever the implementation of any contact spring changes. All of the contact springs must satisfy the following constraints: no tension appears in an open contact and no penetration occurs in a close one. These two conditions can be fulfilled by an openclose iteration process to control the installation and uninstallation of every contact spring in each time step. Readers are recommended to refer to the dissertation [1] for the implementation of openclose iterations.
As illustrated in Figure 6(a), the vertex of the element and the edge of the element are in a contact, where a normal contact spring with the stiffness and a shear contact spring with the stiffness are applied. As shown in Figure 6(b), assume the point on the edge is the contact point. is a vertex of the block which is different from the block where , , and locate. Denote the coordinate of at by , and . The displacement of in the th time step can be expressed asSubstitute (1) into the above equation, and we havewhere is the abbreviation of , .
(a)
(b)
In every iteration step, the normal contact spring and the shear contact spring contribute the strain energies and , respectively, and they can be expressed aswhere and denote the normal and the shear contact displacement, respectively. As shown in Figure 6(b), is the distance from to the edge and can be computed as follows:where is the area of the triangle and denotes the length of the edge . Since is the projection of the line on the edge or its extension line, we haveAs follows, the minimizations of the potential energy items and are derived, respectively.
3.7.1. Normal Contact Submatrices
To obtain the derivative of in (51a), we compute the normal contact displacement firstly. in (52) can be deduced by using (48) asFor a small enough time step, the last item in (54) will be a secondorder infinitesimal and hence can be ignored. Then, we haveIn the same case, can be approximated as :Substituting (50), (55), and (56) into (52) leads toSubstitute (57) into (51a), and then minimize with respect to and :whereThe submatrices in (58a), (58b), (58c), and (58d) are added to the global stiffness matrix in (5), while the submatrices in (58e) and (58f) are added to the right item of the equilibrium equations in (5).
3.7.2. Shear Contact Submatrices
In order to obtain the minimization of the strain energy in (51b), we compute the detailed formulation of the shear contact displacement first. To this end, substitute (48) into (53), and we haveWith a small enough stepsize, the last two items in the equation above are infinitesimals and hence can be ignored, while can be approximately replaced by . Therefore, substituting (50) into (60) leads towhereNow, substitute (61) into (51b) and then minimize with respect to and , and we obtainwhich are added to the global stiffness matrix and the right item in (5), respectively, and here,
3.8. Friction Submatrices
When the MohrCoulomb failure criterion is satisfied along the interface of two contact elements, relative sliding will occur. In this case, a normal contact spring needs to be applied at the contact point to resist the penetration along the sliding surface, and the friction effect has to be taken into account if the friction angle of the joint is not zero. In the BBM, frictions are considered as a kind of external forces acting on the elements on the two sides of the sliding surfaces.
As shown in Figure 7(a), the vertex of the element is sliding along the edge of the element , . A normal contact spring with the stiffness is implemented to attempt to push the element out of the element . As shown in Figure 7(b), the point on is the contact point, and are the sliding displacements of and along the directed edge , respectively, and is the normal contact distance which has the computing formulation given in (57). If the friction angle is not zero, the normal contact force and the friction force can be computed as follows:where the sign function sgn() takes the values of 1, 0, and −1 in the cases of , , and , respectively. As the friction force acts on both sides of the sliding surfaces, its potential energy can be calculated as the following summation:At the time node , , assume the coordinate of the point is , . In this way, the displacements and during the th time step can be computed bywhere denotes the length of at . For a small enough time interval , can be approximated by , which can be computed by (56). In this case, substituting (1) and (3) into (67a) and (67b) leads toSubstitute (68a) and (68b) into (66) and then minimize with respect to and , respectively, and we haveThe above submatrices are added to the right item of (5).
(a)
(b)
4. Numerical Examples
The computer program for the BBM has been coded in the C language, whose flowchart is illustrated in Figure 8. The simultaneous equilibrium equations formed in each time step are solved by the symmetric successive overrelaxation preconditioned conjugate gradient (SSORPCG) method [35]. In this section, three elaborate benchmark numerical examples are conducted to verify the formulations derived in the previous sections. The following stress contour figures are plotted based on nodal stresses, and the stress on every mesh node is the mean stress of the elements which overlap with each other at the grid point.
4.1. Parameter Settings
As shown in Table 1, the parameters used in the experiments are given, including time step interval , total time steps , the error tolerance for time iterations, the SSOR relaxation factor , the error tolerance for SSORPCG iterations, the coefficient used in (4) and (29), the normal contact spring stiffness , and the shear contact spring stiffness .
 
denotes Young’s modulus of the material in Section 4.4. 
4.2. Experiment I
The first experiment is designed to examine the capability of the BBM to analyze elasticity problems. As shown in Figure 9(a), a cantilever beam subjected to a concentrated load KN at the free end is considered. The length and the height of the beam are m and m, respectively. The material parameters of the beam are as follows: Young’s modulus MPa, Poisson’s ratio , and the mass density Kg/m^{2}. The displacement of this problem has the following analytical solution [36]: where is the moment of inertia.
(a)
(b)
The DDA, the BBM, and the PLANE2 element (PLANE2 is defined by six nodes having two degrees of freedom at each node and possesses a quadratic displacement behavior) of the commercial software ANSYS are implemented to solve this problem. The computation model of the DDA contains only one block, that is, the beam, so there are 6 degrees of freedom (DOFs) involved. Both of the computations of the BBM and the FEM adopt the unstructured triangular mesh as shown in Figure 9(b), where there exist totally 2202 triangle elements. The DOFs consumed by the three methods are compared in Table 2, in which the DOFs used by the FEM are two times those by the BBM.

This statics problem is solved through the NewtonRaphson incremental method. Here, time step indicates the iteration counter. The stopping criterion is to exceed the maximum iteration number or to satisfy the following inequality constraint:where time step ; element number ; = denotes the computed displacement at in the th time step; is the square norm; is the error tolerance for NewtonRaphson iterations. In this example, the computation of the BBM converges at .
Comparative Task 1. Figure 10 plots the computed displacements on the nodes distributed along the line = . According to the figure, both of the displacement curves of the BBM and the FEM coincide with that of the analytical solution, while the results obtained by the DDA are low in accuracy. So the accuracy of the original DDA indeed can be significantly improved by our method.
(a)
(b)
Comparative Task 2. In Table 3, the displacement components and at the point computed by the BBM and the FEM are compared with the theoretical solutions. From the table, it can be found that both the BBM and the FEM are very close to the analytical displacements.

Comparative Task 3. Figures 11 and 12 plot the computed contours of the displacement components and throughout the considered cantilever, respectively. According to these pictures, we can find that BBM results agree with ANSYS solutions very well in spite of the introduction of the bonding springs to glue elements.
(a)
(b)
(a)
(b)
Comparative Task 4. As shown in Figure 13, the vonMises stress distributions figured out by the BBM and the FEM with PLANE2 element are illustrated. One can observe that the results obtained by the two methods are very close. Based on this comparative, we can see that the BBM is capable of giving accurate enough stress analysis for this elasticity problem.
(a)
(b)
4.3. Experiment II
The objective of the second example is to demonstrate the capability of the BBM to analyze the elastostatics problem defined on an irregularly shaped domain. As shown in Figure 14(a), a thin plate with a circular hole at the center is pressured by a uniform horizontal stress MPa. Assume that Young’s modulus, Poisson’s ratio, the thickness, and the mass density of the plate are 2900 GPa, 0.3, 0.2 m, and 0.01 Kg/m^{2}, respectively. The radium of the hole is ; the length and the width of the plate are and , respectively, where m, m, and m. Set the origin of the Cartesian coordinate system at the point , the center of the hole, and the axis toward the right. As shown in Figure 14(a), we also use the polar coordinate system , where the pole locates at the origin , is the polar angle in radian, and is the distance in meter away from . In the case of infinite and , this problem has the following analytical solution derived in [37, 38]:
(a)
(b)
We implement both the BBM and PLANE2 element of ANSYS to numerically solve this issue. Onequarter of the domain, that is, the area , is analyzed due to the symmetry, and the symmetrical boundary constraints are applied along the and , respectively. The mesh adopted in both the BBM and the FEM is plotted in Figure 14(b), where 3253 triangles are generated. During the computation, the BBM and the FEM consume totally 19518 DOFs and 39036 DOFs, respectively. Since this example involves only small elastic deformations, it is not necessary to use incremental iterations for the solution. So we use a large time interval and only one time step in the BBM code.
Comparative Task 1. In Table 4, the horizontal stresses computed by the BBM and the FEM are compared with the analytical solutions at the points , , , , and . According to this table, both the BBM and the FEM can give accurate enough horizontal stresses at , , , and , while computed by the two methods have some difference from the theoretical solution. Actually, such errors are caused due to the fact that the plate sizes and are assumed to be infinite in the analytical formulations in (72a) and (72b), but they have the specific values of and in the computational models of the BBM and the FEM. So it still can be noticed that the stress analysis by the BBM agrees with that by the FEM very well.
