Discontinuous Deformation Analysis Enriched by the Bonding Block Model

Sun, Yue; Chen, Qian; Feng, Xiangchu; Wang, Ying

doi:https://doi.org/10.1155/2015/723263

Mathematical Problems in Engineering

On this page

Abstract Introduction Conclusion Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2015 | Article ID 723263 | https://doi.org/10.1155/2015/723263

Discontinuous Deformation Analysis Enriched by the Bonding Block Model

Yue Sun,¹Qian Chen,²Xiangchu Feng,¹and Ying Wang³

Academic Editor: Fabio Tramontana

Received21 Jan 2015

Accepted05 Mar 2015

Published07 Apr 2015

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 open-close 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 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–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 open-close 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–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 nodal-based 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 open-close 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 ill-conditioned 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: vertex-to-vertex, edge-to-edge, and vertex-to-edge contacts. In two-dimensional problems, both vertex-to-vertex and edge-to-edge contacts can be decomposed into couples of vertex-to-edge 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 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 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 second-order 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 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 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 (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.

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 SSOR-PCG iterations, the coefficient used in (4) and (29), the normal contact spring stiffness , and the shear contact spring stiffness .

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². 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 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 ; element number ; = denotes the computed displacement at in the th time step; is the square norm; is the error tolerance for Newton-Raphson 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 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.

(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², 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. One-quarter 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.

Comparative Task 2. In Table 5, the displacement components and computed by the BBM and the FEM are compared at the five measurement points: , , , , and . From the table, we can find that the 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. Since such small displacements (less than 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.

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.

(a)

(b)

(a)

(b)

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.

(a)

(b)

(a)

(b)

4.4. 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).

(a)

(b)

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 s. Below, the computed displacement at the point as marked in Figure 24(a) is concerned. The analytical solution can be deduced based on kinematics theories; that is, and . 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.

5. 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors sincerely thank the referees and the editor for their comments and suggestions. This research is supported by the National Natural Science Foundation of China (Grants nos. 60902098 and 1110222).

References

G.-H. Shi, Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block systems [Ph.D. thesis], University of California, Berkeley, Calif, USA, 1988.
G.-H. Shi and R. E. Goodman, “Generalization of two-dimensional discontinuous deformation analysis for forward modelling,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 13, no. 4, pp. 359–380, 1989.
View at: Publisher Site | Google Scholar
C. T. Lin, B. Amadei, J. Jung, and J. Dwyer, “Extensions of discontinuous deformation analysis for jointed rock masses,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 33, no. 7, pp. 671–694, 1996.
View at: Publisher Site | Google Scholar
S. A. R. Beyabanaki, R. G. Mikola, S. O. R. Biabanaki, and S. Mohammadi, “New point-to-face contact algorithm for 3-D contact problems using the augmented Lagrangian method in 3-D DDA,” Geomechanics and Geoengineering, vol. 4, no. 3, pp. 221–236, 2009.
View at: Publisher Site | Google Scholar
S. A. R. Beyabanaki, B. Ferdosi, and S. Mohammadi, “Validation of dynamic block displacement analysis and modification of edge-to-edge contact constraints in 3-D DDA,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 7, pp. 1223–1234, 2009.
View at: Publisher Site | Google Scholar
Y. Ning, J. Yang, G. Ma, and P. Chen, “Contact algorithm modification of DDA and its verification,” in Analysis of Discontinuous Deformation: New Developments and Applications, G. Ma and Y. Zhou, Eds., pp. 73–81, 2010.
View at: Google Scholar
S. Amir Reza Beyabanaki and A. C. Bagtzoglou, “Non-rigid disk-based DDA with a new contact model,” Computers and Geotechnics, vol. 49, pp. 25–35, 2013.
View at: Publisher Site | Google Scholar
H. Bao, Z. Zhao, and Q. Tian, “On the implementation of augmented lagrangian method in the two-dimensional discontinuous deformation analysis,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 6, pp. 551–571, 2014.
View at: Publisher Site | Google Scholar
J.-H. Wu, “Seismic landslide simulations in discontinuous deformation analysis,” Computers and Geotechnics, vol. 37, no. 5, pp. 594–601, 2010.
View at: Publisher Site | Google Scholar
J.-H. Wu and C.-H. Chen, “Application of DDA to simulate characteristics of the tsaoling landslide,” Computers and Geotechnics, vol. 38, no. 5, pp. 741–750, 2011.
View at: Publisher Site | Google Scholar
Y. Zhang, X. Fu, and Q. Sheng, “Modification of the discontinuous deformation analysis method and its application to seismic response analysis of large underground caverns,” Tunnelling and Underground Space Technology, vol. 40, pp. 241–250, 2014.
View at: Publisher Site | Google Scholar
C. J. Pearce, A. Thavalingam, Z. Liao, and N. Bićanić, “Computational aspects of the discontinuous deformation analysis framework for modelling concrete fracture,” Engineering Fracture Mechanics, vol. 65, no. 2-3, pp. 283–298, 2000.
View at: Publisher Site | Google Scholar
Y.-J. Ning, J. Yang, X. M. An, and G. W. Ma, “Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework,” Computers and Geotechnics, vol. 38, no. 1, pp. 40–49, 2011.
View at: Publisher Site | Google Scholar
Q. Tian, Z. Zhao, and H. Bao, “Block fracturing analysis using nodal-based discontinuous deformation analysis with the double minimization procedure,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 9, pp. 881–902, 2013.
View at: Publisher Site | Google Scholar
Y. X. Ben, Y. Wang, and G.-H. Shi, “Development of a model for simulating hydraulic fracturing with DDA,” in Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation (ICADD '13), pp. 169–175, Fukuoka, Japan, August 2013.
View at: Google Scholar
Y.-Y. Jiao, H.-Q. Zhang, X.-L. Zhang, H.-B. Li, and Q.-H. Jiang, “A two-dimensional coupled hydromechanical discontinuum model for simulating rock hydraulic fracturing,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 39, no. 5, pp. 457–481, 2015.
View at: Publisher Site | Google Scholar
C. Y. Koo and J. C. Chern, “The development of DDA with third order displacement function,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, M. R. Salami and D. Banks, Eds., pp. 342–349, TSI Press, Berkeley, Calif, USA, 1996.
View at: Google Scholar
M. Y. Ma, M. Zaman, and J. H. Zhu, “Discontinuous deformation analysis using the third order displacement function,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, M. R. Salami and D. Banks, Eds., pp. 383–394, TSI Press, Berkeley, Calif, USA, June 1996.
View at: Google Scholar
S. M. Hsiung, “Discontinuous deformation analysis (DDA) with nth order polynomial displacement functions,” in Proceedings of the 38th U.S. Rock Mechanics Symposium on Rock Mechanics in the National Interest, pp. 1413–1420, American Rock Mechanics Association, Balkema, Rotterdam, The Netherlands, Washington, DC, USA, 2001.
View at: Google Scholar
S. A. R. Beyabanaki, A. Jafari, and M. R. Yeung, “High-order three-dimensional discontinuous deformation analysis (3-D DDA),” International Journal for Numerical Methods in Biomedical Engineering, vol. 26, no. 12, pp. 1522–1547, 2010.
View at: Publisher Site | Google Scholar
A. Q. Wu, Y. Zhang, and S.-Z. Lin, “Complete and high order polynomial displacement approximation and its application to elastic mechanics analysis based on DDA,” in Advances in Discontinuous Numerical Methods and Applications in Geomechanics and Geoengineering, T. Sasaki, Ed., pp. 43–54, 2012.
View at: Google Scholar
B. Lu, A. Wu, and X. Ding, “Mixed higher-order discontinuous deformation analysis,” in Frontiers of Discontinuous Numerical Methods and Practical Simulations in Engineering and Disaster Prevention, G. Chen, Y. Ohnishi, L. Zheng, and T. Sasaki, Eds., pp. 243–248, CRC Press, 2013.
View at: Google Scholar
T. C. Ke, “Artificial joint based DDA,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, M. R. Salami and D. Banks, Eds., pp. 326–333, TSI Press, Albuquerque, NM, USA, 1996.
View at: Google Scholar
Y. H. Zhang and Y. M. Cheng, “Coupling of FEM and DDA methods,” The International Journal of Geomechanics, vol. 2, no. 4, pp. 503–517, 2002.
View at: Publisher Site | Google Scholar
M. S. Khan, Investigation of discontinuous deformation analysis for application in jointed rock masses [Ph.D. thesis], Department of Civil Engineering of University of Toronto, 2010.
Z.-Y. Zhao and J. Gu, “Stress recovery procedure for discontinuous deformation analysis,” Advances in Engineering Software, vol. 40, no. 1, pp. 52–57, 2009.
View at: Publisher Site | Google Scholar
M. Zhang, H. Yang, and Z. K. Li, “A coupling model of the discontinuous deformation analysis method and the finite element method,” Tsinghua Science and Technology, vol. 10, no. 2, pp. 221–226, 2005.
View at: Publisher Site | Google Scholar
K. K. Shyu, Nodal-based discontinuous deformation analysis [Ph.D. thesis], Department of Civil Engineering, University of California, Berkeley, Calif, USA, 1993.
Q. Chang, Non-linear dynamic discontinuous deformation analysis with finite element meshed block systems [Ph.D. thesis], University of California, Berkeley, Calif, USA, 1994.
R. Grayeli and A. Mortazavi, “Discontinuous deformation analysis with second-order finite element meshed block,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 30, no. 15, pp. 1545–1561, 2006.
View at: Publisher Site | Google Scholar
R. Grayeli and K. Hatami, “Implementation of the finite element method in the three-dimensional discontinuous deformation analysis (3D-DDA),” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 32, no. 15, pp. 1883–1902, 2008.
View at: Publisher Site | Google Scholar
S. A. R. Beyabanaki, A. Jafari, S. O. R. Biabanaki, and M. R. Yeung, “Nodal-based three-dimensional discontinuous deformation analysis (3-D DDA),” Computers and Geotechnics, vol. 36, no. 3, pp. 359–372, 2009.
View at: Publisher Site | Google Scholar
H. R. Bao and Z. Y. Zhao, “Modeling brittle fracture with the nodal-based discontinuous deformation analysis,” International Journal of Computational Methods, vol. 10, no. 6, Article ID 1350040, 26 pages, 2013.
View at: Publisher Site | Google Scholar | MathSciNet
G.-H. Shi, “Modeling dynamic rock failure by discontinuous deformation analysis with simplex integrations,” in Proceedings of the 1st North American Rock Mechanics Symposium, pp. 591–598, Austin, Tex, USA, 1994.
View at: Google Scholar
X. Chen and K. K. Phoon, “Some numerical experiences on convergence criteria for iterative finite element solvers,” Computers and Geotechnics, vol. 36, no. 8, pp. 1272–1284, 2009.
View at: Publisher Site | Google Scholar
C. E. Augarde and A. J. Deeks, “The use of Timoshenko's exact solution for a cantilever beam in adaptive analysis,” Finite Elements in Analysis and Design, vol. 44, no. 9-10, pp. 595–601, 2008.
View at: Publisher Site | Google Scholar
S. P. Timoshenko and J. N. Goodier, Theory of Elasticity, McGraw-Hill, New York, NY, USA, 3rd edition, 1970.
J. C. Jaeger, N. G. W. Cook, and R. W. Zimmerman, Fundamentals of Rock Mechanics, John Wiley & Sons, 4th edition, 2007.

Copyright

Copyright © 2015 Yue Sun et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1225

Downloads

1011

Citations