Discontinuous Deformation Analysis Enriched by the Bonding Block Model

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 ourmethod, 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.


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 open-close iterations.In the literature [3][4][5][6][7][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][10][11], crack propagations [12][13][14], and hydraulic fractures [15,16].
However, the DDA employs the first-order 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 higher-order polynomial functions to approximate block displacements [17][18][19][20][21][22].But the accuracy of this method is only adequate for regularly shaped blocks.Another way is the so-called 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 so-called 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 nodal-based DDA [28][29][30][31][32][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 open-close 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 open-close iterations along artificial joints.Compared with the coupling DDA-FEM 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.

Concepts of the BBM
2.1.Subdivision.For a jointed rock mass, assume it contains Λ blocks occupying the volumes Ω 1 , Ω 2 , . . ., Ω Λ .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 I  , I +1 , . . ., I  .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 I 1 , I 2 , . . ., I  .On each element, say I  ,  = 1, 2, . . ., , 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 I  ; {  } is composed of the unknowns to be determined of the element; that is, where  0 , V 0 are the translational displacements of the barycenter,  0 is the rotation of the centroid, and   ,   , and   denote the strain components of I  .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.

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, I  and I  are two adjacent elements of the block Ω  .From Figure 3(a), it can be recognized that  2  3 and  4  5 are boundary edges, while  1  2 ,  3  1 ,  5  6 , and  6  4 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 open-close iterations.
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  3  1 and  6  4 are glued together by two bonding springs applied between  1 and  6 and between  3 and  4 , 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 I  and I  : where κ is a constant,  is the average shear modulus of the two elements, and   and   denote the areas of I  and I  , 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  ≈ 1) to 1000 (for  ≫ 1).
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 ill-conditioned global stiffness matrix.For the simplicity, we use  instead of   to denote bonding spring stiffness in the following sections.

Computing Formulations
The BBM implements an incremental iteration process to solve both dynamics and statics problems of the element system {I 1 , I 2 , . . ., I  } on a given time mesh {0 =  0 <  1 < ⋅ ⋅ ⋅ <   < ⋅ ⋅ ⋅ <   }.In particular, time steps used in statics problems represent the loading steps.Below, the th time step indicates the interval [ −1 ,   ].Obviously, it is equivalent to say the value of a variable at the end of the ( − 1)th step and the one at the beginning of the th step.
For the element system {I 1 , I 2 , . . ., I  } 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 6 × 6 and 6 × 1, respectively, while {   } denotes the unknown vector to be determined, ,  = 1, 2, . . ., .The detailed formulations of these submatrices, [  ] and {  }, are derived as follows.
3.1.Element Elastic Submatrices.The strain energy of the th element I  can be expressed as where 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 as Substituting ( 7) and ( 8) into (6), we have Minimizing Π   with respect to {   } leads to which are added to the stiffness matrix and the right item of (5), respectively.

Element Inertia Submatrices.
The potential energy contributed by the inertia force acting on the element I  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 I  .Substituting (1) into (11) leads to Using the Newmark time integrating method, we have the following finite difference approximation: where Δ  =   −  −1 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 have Minimizing the right items of the equation above with respect to 2 Δ  (∬ where 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.

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  1 and  2 are a pair of initially overlapped vertexes of the elements I  and I  which are within the same block.A bonding spring with the stiffness  is applied between  1 and  2 to glue I  and I  together.Denote the coordinate of  1 and  2 at the time node   by (  1 ,   1 ) and (  2 ,   2 ), respectively.Firstly, we compute the elongation   of the bonding spring between  1 and  2 at the time node   .As illustrated in Figure 4(b),   is actually the length of   →  1  2 and can be computed from the following equality: where {  } denotes the vector   →  1  2 at   , that is, Assume the displacements of  1 and  2 in the th time step are (  1 , V  1 ) and (  2 , V  2 ), respectively.Then, the coordinates of  1 and  2 at   can be decomposed into the incremental form below: Substituting (20a) and (20b) into ( 19) leads to Substitute ( 3) and (1) into the equation above, and then we have where 22) into (18), then  2  can be expanded as the following summation: Next, we compute the strain energy Π   of the bonding spring between  1 and  2 in the th time step [ −1 ,   ].Π   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 { −1 } is known in the th time step, the derivative of  0 is 0.
Minimize   ,   ,   , and   with respect to {  } and {  }, and we have These 6 × 6 submatrices in the equation above are added into the global stiffness matrix in (5).Minimizing   and   with respect to {  } and {  } leads to which are added to the right side of (5).

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 I  is a fixed point,  = 1, 2, . . ., .Denote the coordinates of  at the time node   by (   ,    ),  = 1, 2, . . ., .Assume there exists a fictitious element I 0 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  0 of I 0 overlaps with the fixed point  at the initial time  0 .In order to resist the displacement of the point , a bonding spring with the stiffness  0 is applied between  and  0 .The stiffness  0 can be determined by the following formulation: where κ is the same factor used in (4).As shown in Figure 5(b), the bonding spring between  and  0 has the elongation   at the time node   .Denote the coordinate of  0 by (, ), which keeps constant during a computation.Since   is actually the length of   →  0 , we have where {   } denotes the vector   →  0  at   ; that is, Assume the displacement of  in the th time step is (   , V   ).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 have where Then we can express the strain energy of the bonding spring between  and  0 as where Obviously, the minimization of  0 with respect to {   } is zero.Then, we minimize   and   : where [  ] and {  } are added to the stiffness matrix and the right item of (5).

Line Loading Submatrices.
Assume a loading {} = (  (, ),   (, )) T is applied along a line upon the boundary of the element I  ,  = 1, 2, . . ., .The coordinate of any point on the line at the time moment   can be determined by the following parametric equations: where (  1 ,   1 ) and (  2 ,   2 ) are the coordinates of the ending points of the line, whose length at the end of the last time step can be computed by where time step  = 1, 2, . . ., .Then, we can obtain the potential energy contributed by the loading {}: V  ( () ,  ()) ) Substituting ( 3) and ( 1) into (40) leads to Minimize Π  with respect to {   }, and we have The 6 × 1 submatrix is added to the right item of the global equation (5).If {} = (  ,   ), which is constant, (42) has the analytical form: where x =  2 +  1 , ỹ =  2 +  1 and ( 0 ,  0 ) is the coordinate of the barycenter of the element I  .
3.6.Volume Force Submatrix.When the element I  bears the body loading (, ) = (  (, ),   (, )) T ,  = 1, 2, . . ., , the potential energy Π  due to the force  is calculated as Substituting (1) into the equations above leads to Minimize Π  with respect to {   }, and we have The 6 × 1 submatrix above is added to the right side of the global simultaneous equation (5).When the body force is constant, say (, ) = (  ,   ) T , (46) becomes where   denotes the area of the element I  .

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: vertexto-vertex, edge-to-edge, and vertex-to-edge contacts.In twodimensional problems, both vertex-to-vertex and edge-toedge contacts can be decomposed into couples of vertex-toedge ones.The BBM adopts the "rigid" contact model.In order to prevent penetration and relative sliding along the (a) 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 open-close 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 open-close iterations.
As illustrated in Figure 6(a), the vertex  1 of the element I  and the edge  2  3 of the element I  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  0 on the edge  2  3 is the contact point. 1 is a vertex of the block which is different from the block where  0 ,  2 , and  3 locate.Denote the coordinate of   at   by (   ,    ),  = 0, 1, 2, 3 and  = 1, 2, . . ., .The displacement of   in the th time step can be expressed as Substitute (1) into the above equation, and we have where [   ] is the abbreviation of [  (  ,   )],  = 0, 1, 2, 3.In every iteration step, the normal contact spring and the shear contact spring contribute the strain energies Π ⊥ and Π ‖ , respectively, and they can be expressed as where ℎ ⊥ and ℎ ‖ denote the normal and the shear contact displacement, respectively.As shown in Figure 6(b), ℎ ⊥ is the distance from  1 to the edge  2  3 and can be computed as follows: where   is the area of the triangle Δ  1  2  3 and   denotes the length of the edge  2  3 .Since ℎ ‖ is the projection of the line  0  1 on the edge  2  3 or its extension line, we have As follows, the minimizations of the potential energy items Π ⊥ and Π ‖ are derived, respectively.

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 have With 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  −1 .Therefore, substituting (50) into (60) leads to where ) . (62) Now, substitute (61) into (51b) and then minimize Π ‖ with respect to {   } and {   }, and we obtain which are added to the global stiffness matrix and the right item in (5), respectively, and here, 3.8.Friction Submatrices.When the Mohr-Coulomb 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  1 of the element I  is sliding along the edge  2  3 of the element I  , ,  = 1, 2, . . ., .A normal contact spring with the stiffness  ⊥ is implemented to attempt to push the element I  out of the element I  .As shown in Figure 7(b), the point  0 on  2  3 is the contact point,  0 and  1 are the sliding displacements of  0 and  1 along the directed edge   →  2  3 , 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 N and the friction force F can be computed as follows: where the sign function sgn() takes the values of 1, 0, and −1 in the cases of  > 0,  = 0, and  < 0, respectively.As the friction force F acts on both sides of the sliding surfaces, its potential energy Π F can be calculated as the following summation: At the time node   ,  = 1, 2, . . ., , assume the coordinate of the point   is (   ,    ),  = 0, 1, 2, 3.In this way, the displacements  0 and  1 during the th time step can be computed by ) , (67a) where   denotes the length of  2  3 at   .For a small enough time interval [ −1 ,   ],   can be approximated by  −1 , which can be computed by (56).In this case, substituting ( 1) and ( 3) into (67a) and (67b) leads to ) , (68a) Substitute (68a) and (68b) into (66) and then minimize Π F with respect to {   } and {   }, respectively, and we have The above 6×1 submatrices are added to the right item of (5).

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 (SSOR-PCG) 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.interval Δ, total time steps , the error tolerance   for time iterations, the SSOR relaxation factor , the error tolerance  for SSOR-PCG iterations, the coefficient κ used in ( 4) and ( 29), the normal contact spring stiffness  ⊥ , and the shear contact spring stiffness  ‖ .

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  = 1 KN at the free end is considered.The length and the height of the beam are  = 10 m and ℎ = 1 m, respectively.The material parameters of the beam are as follows: Young's modulus  = 300 MPa, Poisson's ratio ] = 0.3, and the mass density  = 0.01 Kg/m 2 .The displacement of this problem has the following analytical solution [36]: where  = ℎ 3 /12 is the moment of inertia.
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 Newton-Raphson 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  = 1, 2, . . ., 100; element number  = 2202; {  } = {  (, )} denotes the computed displacement at (, ) in the th time step; ‖ ⋅ ‖ 2 is the square norm;   = 10 −8 is the error tolerance for Newton-Raphson iterations.In this example, the computation of the BBM converges at  = 3.
Comparative Task 1. Figure 10 plots the computed displacements on the nodes distributed along the line  = {(, ) : 0 ⩽  ⩽ ,  = −ℎ/2}.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.
Comparative Task 2. In Table 3, the displacement components   and   at the point (, −ℎ/2) 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.x (m) Comparative Task 4. As shown in Figure 13, the von-Mises 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.

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  ∞ = 1 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 2; the length and the width of the plate are 2 and 2, respectively, where  = 0.5 m,  = 10 m, and  = 5 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]: We implement both the BBM and PLANE2 element of ANSYS to numerically solve this issue.One-quarter of the Table 4: The stresses   computed by the BBM and the FEM are compared with the analytical solutions at the points (0, 0.5), (0.5, 0), (10, 0),  (10,5), and (0, 5).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.

Stress/MPa
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 (0, 0.5), (0.5, 0), (10, 0),  (10,5), and (0, 5).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  = 10 and  = 5 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.Comparative Task 2. In Table 5, the displacement components   and   computed by the BBM and the FEM are compared at the five measurement points: (0, 0.5), (0.5, 0), (10, 0),  (10,5), and (0, 5).From the table, we can find that the The FEM result The BBM result displacements solved out by the BBM are very close to the FEM results.The subtle differences between the results of the two methods on   (),   (),   (), and   () are due to the fact that the BBM imports the symmetric boundary constraints by using penalty-like bonding springs.The FEM result

Analytical solution
The BBM result   Since such small displacements (less than 10 −11 m) can be interpreted physically as the local elastoplastic deformation of the plate near the partially fixed boundaries, the relative errors between our method and the FEM on the computation of   (),   (),   (), and   () are acceptable.Comparative Task 3.Then, we compare the computed deformations along the and -axis with the theoretical results.Figure 15 plots the displacements   and   obtained from the BBM and the FEM results on the nodes along the boundary .And, in Figure 16,   and   solved by the two methods on the nodes along  are drawn.From the figures, it can be observed that the BBM results are in accord with the FEM solutions.

Mathematical Problems in Engineering
Comparative Task 4. As shown in Figures 17, 18, and 19, the computed horizontal stresses   on the nodes along , , and the arc  are illustrated, respectively.According to the three plots, one can find that both of the results of the BBM and the FEM are consistent with the analytical solutions.
Comparative Task 5. Next, we compare the BBM and the FEM on analyzing the deformation distributions of the quarter of the plate.As shown in Figures 20 and 21, the computed contour figures of the stresses   throughout the domain are illustrated, and it can be observed that the BBM results agree with the FEM very well.
Comparative Task 6.In Figures 22 and 23, the distributions of the displacements obtained by the BBM and the FEM are demonstrated.According to these contour plots, we can see that the result of the BBM conforms with the solution of the FEM.

Experiment III.
The third numerical example is to demonstrate the accuracy and the effectiveness of the BBM in the analysis of a frictionless sliding problem.As shown in Figure 24(a), a small rectangular block subjected to a horizontal traction force  slides along the top surface of a fixed flat plate.The origin of the coordinate system locates at the plate centroid, and the  direction points to the right.Both the upper block and the lower plate are all linear isotropic elastic and are with the material constants in Figure 24(a).We implement the BBM to solve this problem by using the quadrilateral mesh as plotted in Figure 24(b), where the upper rectangular block and the lower flat plate are subdivided into 8 and 160 elements, respectively.In the computation, set the uniform time step length as 10 −6 s.Below, the computed displacement at the point (−1, 0.3) as marked in Figure 24

Conclusion
In this paper, we extended the concept of DDA blocks, which are separated by sets of geological joints.Our method reconstructs a block system as an assembly of block-like elements through an advanced subdivision process.These elements can be triangles, quadrilaterals, and even the original DDA blocks without partitioning.Overlapped element edges within a block are glued together by bonding springs, while contact springs are implemented on the element edges that compose the boundaries of original blocks.In this way, block displacement and stress fields can be refined, and the efficient contact algorithms developed in the standard DDA are also retained.Such scheme is easy to use and can be adopted in both the 2-dimensional and the 3-dimensional DDA although it was performed in the 2-dimensional case only here.Numerical experiments indicate that the proposed method is accurate and effective for both continuous and discontinuous problems.Therefore, this procedure is potential to be applied in crack propagation analysis.However, such problem exceeds the main concern of this work.We will continue the task in the foreseeable future.

Figure 1 :
Figure 1: Block subdivision adopted in the BBM.

Figure 3 :
Figure 3: Overlapped element edges within a block are glued together by a pair of bonding springs applied between the coincided vertexes.(a) The initial status of bonding springs.(b) Bonding spring elongations.

Figure 4 :
Figure 4: Bonding springs applied along overlapped element edges.(a) Bonding spring length is zero at  0 .(b) The bonding spring between  1 and  2 has the elongation   at   .

MathematicalFigure 5 :
Figure 5: Bonding springs are applied at fixed points.(a) The length of the bonding spring applied between the fixed point  of the element I  and the vertex  0 of the fictitious support element I 0 is zero at  0 .(b) At   , the bonding spring between  0 and  has the elongation   .

Figure 6 :
Figure 6: Contact springs between the elements I  and I  .(a) The normal contact spring with the stiffness  ⊥ and the shear contact spring with the stiffness  ‖ .(b) The normal contact distance ℎ ⊥ and the shear contact distance ℎ ‖ .

Figure 7 :Figure 8 :
Figure 7: The vertex  1 of the element I  is sliding along the edge  2  3 of the element I  .(a) A normal contact spring with the stiffness  ⊥ is applied (b).

Figure 9 :
Figure 9: A cantilever subjected to a concentrated load at the free end.(a) Problem illustration on the geometry configuration and the given mechanical parameters.(b) Triangular mesh used in both the BBM and the FEM.

Figure 10 :Figure 11 :Figure 12 :
Figure 10: The computed displacement components   and   of the nodes along the line .

Figure 20 :Figure 21 :
Figure 20: The contour plots of the computed distributions of the horizontal stress   (Pa) throughout the plate.(a) The result of the BBM.(b) The result of the FEM.

Figure 22 :
Figure 22: The contour plots of the computed distributions of the displacement component   (m) throughout the plate.(a) The result of the BBM.(b) The result of the FEM.

Figure 23 :Figure 24 :
Figure 23: The contour plots of the computed distributions of the displacement component   (m) throughout the plate.(a) The result of the BBM.(b) The result of the FEM.
(a) is concerned.The analytical solution ( Analyt  ,  Analyt  ) can be deduced based on kinematics theories; that is,  Analyt  = 50000 ⋅  2 and  Analyt  = 0.As shown in Figure 25, the computed displacement components   and   versus the time are plotted.From the figure, it can be observed that the results of the BBM agree with the analytical solutions very well.

Figure 25 :
Figure 25: The computed displacements measured at the point  in the  and  directions, respectively.

Figure 2 :
The subdivision process implemented on the block Ω  .

Table 1
, the parameters used in the experiments are given, including time step

Table 1 :
The values of the parameters used in the experiments. † 10 2  †  denotes Young's modulus of the material in Section 4.4.

Table 2 :
Comparing the DOFs involved in the DDA, the FEM, and the BBM for Section 4.2.

Table 3 :
The displacements at (, −ℎ/2) obtained by the analytical solution, the BBM, and the FEM.

Table 5 :
The displacements   and   computed by the BBM and