# Discontinuous Deformation Analysis Enriched by the Bonding Block Model

Authors: Yue Sun; Qian Chen; Xiangchu Feng; Ying Wang
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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 Figure1, 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.Figure 1
Block subdivision adopted in the BBM.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 Section3, 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 Ω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 Ii, Ii+1,…,Ij. As illustrated in Figure 2(c), the overlapped element edges are separated from each other. Proceed such subdivision over all blocks, and assume M elements are generated, including I1,I2,…,IM. On each element, say Ii,i=1,2,…,M, the displacement at any point (x,y) can be approximated by the following function:(1)Ux,y=Tix,yDi,where [Ti(x,y)] is the displacement mode matrix:(2)Tix,y=10y¯-yx-x¯0y-y¯201x-x¯0y-y¯x-x¯2,in which (x¯,y¯) denotes the barycenter coordinate of the element Ii; {Di} is composed of the unknowns to be determined of the element; that is,(3){Di}=u0v0r0εxεyγxyT,where u0, v0 are the translational displacements of the barycenter, r0 is the rotation of the centroid, and εx, εy, and γxy denote the strain components of Ii.Figure 2
The subdivision process implemented on the blockΩα.
(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 Section3 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 Figure2, Ii and Ij are two adjacent elements of the block Ωα. From Figure 3(a), it can be recognized that P2P3¯ and P4P5¯ are boundary edges, while P1P2¯, P3P1¯, P5P6¯, and P6P4¯ 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.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.
(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 Figure3(b), the overlapped edges P3P1¯ and P6P4¯ are glued together by two bonding springs applied between P1 and P6 and between P3 and P4, 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 Ii and Ij:(4)κij=κ~G1/2si+sj,where κ~ is a constant, G is the average shear modulus of the two elements, and si and sj denote the areas of Ii and Ij, 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 κij to denote bonding spring stiffness in the following sections.


## 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 Ii, Ii+1,…,Ij. As illustrated in Figure 2(c), the overlapped element edges are separated from each other. Proceed such subdivision over all blocks, and assume M elements are generated, including I1,I2,…,IM. On each element, say Ii,i=1,2,…,M, the displacement at any point (x,y) can be approximated by the following function:(1)Ux,y=Tix,yDi,where [Ti(x,y)] is the displacement mode matrix:(2)Tix,y=10y¯-yx-x¯0y-y¯201x-x¯0y-y¯x-x¯2,in which (x¯,y¯) denotes the barycenter coordinate of the element Ii; {Di} is composed of the unknowns to be determined of the element; that is,(3){Di}=u0v0r0εxεyγxyT,where u0, v0 are the translational displacements of the barycenter, r0 is the rotation of the centroid, and εx, εy, and γxy denote the strain components of Ii.Figure 2
The subdivision process implemented on the blockΩα.
(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 Section3 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 Figure2, Ii and Ij are two adjacent elements of the block Ωα. From Figure 3(a), it can be recognized that P2P3¯ and P4P5¯ are boundary edges, while P1P2¯, P3P1¯, P5P6¯, and P6P4¯ 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.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.
(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 Figure3(b), the overlapped edges P3P1¯ and P6P4¯ are glued together by two bonding springs applied between P1 and P6 and between P3 and P4, 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 Ii and Ij:(4)κij=κ~G1/2si+sj,where κ~ is a constant, G is the average shear modulus of the two elements, and si and sj denote the areas of Ii and Ij, 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 κij 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{I1,I2,…,IM} on a given time mesh {0=t0<t1<⋯<tn<⋯<tN}. In particular, time steps used in statics problems represent the loading steps. Below, the nth time step indicates the interval [tn-1,tn]. Obviously, it is equivalent to say the value of a variable at the end of the (n-1)th step and the one at the beginning of the nth step.For the element system{I1,I2,…,IM} 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 nth step leads to the global simultaneous equilibrium equations:(5)K11K12⋯K1MK21K22⋯K2M⋮⋮⋱⋮KM1KM2⋯KMMD1nD2n⋮DMn=F1F2⋮FM,where the submatrices [Kij] and {Fi} are of the sizes 6×6 and 6×1, respectively, while {Din} denotes the unknown vector to be determined, i,j=1,2,…,M. The detailed formulations of these submatrices, [Kij] and {Fi}, are derived as follows.


### 3.1. Element Elastic Submatrices

The strain energy of theith element Ii can be expressed as(6)Πei=12∬Ii‍εinTσin+2σ¯indxdy,where the vectors {εin}, {σin}, and {σ¯in} denote the strain, stress, and the initial stress in the nth time step, respectively. For isotropic linearly elastic materials, {σin} and {εin} have the following relationship:(7)σin=E1-ν21ν0ν10001-ν2εin=Eiεin,in which E and ν indicate Young’s modulus and Poisson’s ratio, respectively. The strain {εin} can be expressed by the displacement as(8)εin=000100000010000001Din=ΘDin.Substituting (7) and (8) into (6), we have(9)Πei=12DinT∬Ii‍ΘTEiΘDin+DinT∬Ii‍σ¯indxdy.Minimizing Πei with respect to {Din} leads to(10)∂2Πei∂DinT∂Din=∬Ii‍ΘTEiΘdxdy⟶Kii,-∂Πei0∂DinT=-∬Ii‍ΘTσ¯indxdy⟶Fi,which 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 elementIi at the moment t can be expressed as follows:(11)Πρi=∬Ii‍Uit,x,yTρ∂2∂t2Uit,x,ydxdy,where {Ui(t,x,y)} 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 Ii. Substituting (1) into (11) leads to(12)Πρi=DitT∬Ii‍ρTiTTi∂2Dit∂t2dxdy.Using the Newmark time integrating method, we have the following finite difference approximation:(13)∂2Dit∂t2=2DinΔn2-2V¯in-1Δn,where Δn=tn-tn-1 is the nth stepsize and V¯in-1 indicates the initial velocity of the nth time step and can be updated by the iteration formulation below:(14)V¯in=2Din-1Δn-V¯in-1,n⩾2.Substituting (13) into (12), we then have(15)Πρi=2DinTΔn2∬Ii‍ρTiTTidxdyDin-2DinTΔn∬Ii‍ρTiTTidxdyV¯in-1.Minimizing the right items of the equation above with respect to {Din} leads to(16)2Δn2∬Ii‍ρTiTTidxdy⟶Kii,(17)2Δn∬Ii‍ρTiTTidxdyV¯in-1⟶Fi,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.


### 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 Figure4(a), the points P1 and P2 are a pair of initially overlapped vertexes of the elements Ii and Ij which are within the same block. A bonding spring with the stiffness κ is applied between P1 and P2 to glue Ii and Ij together. Denote the coordinate of P1 and P2 at the time node tn by (x1n,y1n) and (x2n,y2n), respectively.Figure 4
Bonding springs applied along overlapped element edges. (a) Bonding spring length is zero att0. (b) The bonding spring between P1 and P2 has the elongation dn at tn.
(a)
(b)Firstly, we compute the elongationdn of the bonding spring between P1 and P2 at the time node tn. As illustrated in Figure 4(b), dn is actually the length of P1P2→ and can be computed from the following equality:(18)dn2=lnTln,where {ln} denotes the vector P1P2→ at tn, that is,(19)ln=x2ny2n-x1ny1n.Assume the displacements of P1 and P2 in the nth time step are (u1n,v1n) and (u2n,v2n), respectively. Then, the coordinates of P1 and P2 at tn can be decomposed into the incremental form below:(20a)x1n,y1n=x1n-1+u1n,y1n-1+v1n,(20b)x2n,y2n=x2n-1+u2n,y2n-1+v2n.Substituting (20a) and (20b) into (19) leads to(21)ln=ln-1+u1nv1n-u2nv2n.Substitute (3) and (1) into the equation above, and then we have(22){ln}=ln-1+Ti1Din-Tj2Djn,where [Ti1]=[Ti(x1,y1)] and [Tj2]=[Tj(x2,y2)]. Substituting (22) into (18), then dn2 can be expanded as the following summation:(23)dn2=ln-1Tln-1+DinTTi1TTi1Din+DjnTTj2TTj2Djn-DinTTi1TTj2Djn-DjnTTj2TTi1Din+ln-1TTi1Din+Ti1TDinTln-1-ln-1TTj2Djn-Tj2TDjnTln-1.Next, we compute the strain energyΠgκ of the bonding spring between P1 and P2 in the nth time step [tn-1,tn]. Πgκ can be expressed as follows:(24)Πgκ=12κdn2.Substitute (23) into the expression above, and Πgκ becomes the following sum:(25)Πgκ=π0+πi+πj+πii+πij+πji+πjj,where(26a)π0=κ2ln-1Tln-1,(26b)πii=κ2DinTTi1TTi1Din,(26c)πjj=κ2DjnTTj2TTj2Djn,(26d)πij=-κ2DinTTi1TTj2Djn,(26e)πji=-κ2DjnTTj2TTi1Din,(26f)πi=κ2ln-1TTi1Din+κ2Ti1TDinTln-1,(26g)πj=-κ2ln-1TTj2Djn-κ2Tj2TDjnTln-1.Then, minimizeΠgκ with respect to the unknowns. As {ln-1} is known in the nth time step, the derivative of π0 is 0. Minimize πii, πij, πji, and πjj with respect to {Di} and {Dj}, and we have(27a)∂2πii∂DinT∂Din=κTi1TTi1⟶Kii,(27b)∂2πij∂DinT∂Djn=κTi1TTj2⟶Kij,(27c)∂2πji∂DjnT∂Din=κTj2TTi1⟶Kji,(27d)∂2πjj∂DjnT∂Djn=κTj2TTj2⟶Kjj.These 6×6 submatrices in the equation above are added into the global stiffness matrix in (5). Minimizing πi and πj with respect to {Di} and {Dj} leads to(28a)-∂πi0∂Di=-κTi1Tln-1⟶Fi,(28b)-∂πj({0})∂Dj=κTj2Tln-1⟶Fj,which 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 Figure5(a), the point P of the element Ii is a fixed point, i=1,2,…,M. Denote the coordinates of P at the time node tn by (xPn,yPn), n=1,2,…,N. Assume there exists a fictitious element I0 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 P0 of I0 overlaps with the fixed point P at the initial time t0. In order to resist the displacement of the point P, a bonding spring with the stiffness κ0 is applied between P and P0. The stiffness κ0 can be determined by the following formulation:(29)κ0=κ~E,where κ~ is the same factor used in (4).Figure 5
Bonding springs are applied at fixed points. (a) The length of the bonding spring applied between the fixed pointP of the element Ii and the vertex P0 of the fictitious support element I0 is zero at t0. (b) At tn, the bonding spring between P0 and P has the elongation dn.
(a)
(b)As shown in Figure5(b), the bonding spring between P and P0 has the elongation dn at the time node tn. Denote the coordinate of P0 by (x¯,y¯), which keeps constant during a computation. Since dn is actually the length of P0P→, we have(30)dn2=lPnTlPn,where {lPn} denotes the vector P0P→ at tn; that is,(31)lPn=xPnyPn-x¯y¯.Assume the displacement of P in the nth time step is (uPn,vPn). Then, the coordinates of P at tn can be written as xPn,yPn=xPn-1+uPn,yPn-1+vPn, which is substituted into (31); that is,(32)lPn=xPn-1+uPnyPn-1+vPn-x¯y¯=ln-1+uPnvPn.Substitute (3) and (1) into the equation above, and then we have(33)lPn=lPn-1+TiPDin,where [TiP] indicates [Ti(xPn,yPn)]. Substituting (33) into (30) leads to(34)dn2=lPn-1TlPn-1+DinTTi1TTi1Din+lPn-1TTi1Din+Ti1TDinTlPn-1.Then we can express the strain energy of the bonding spring between P and P0 as(35)ΠPκ0=12κ0dn2=π0+πii+πi,where(36a)π0=κ02lPn-1TlPn-1,(36b)πii=κ02DinTTi1TTi1Din,(36c)πi=κ0Ti1TDinTlPn-1.Obviously, the minimization of π0 with respect to {Din} is zero. Then, we minimize πii and πi:(37a)∂2πii∂DinT∂Din=κ0TPTTP⟶Kii,(37b)∂πi∂DinT=κ0lPn-1TTP⟶Fi,where [Kii] and {Fi} are added to the stiffness matrix and the right item of (5).


### 3.5. Line Loading Submatrices

Assume a loading{F}=(Fx(x,y),Fy(x,y))T is applied along a line upon the boundary of the element Ii, i=1,2,…,M. The coordinate of any point on the line at the time moment tn can be determined by the following parametric equations:(38)xn=xnξ=x2n-x1nξ+x1n,yn=ynξ=y2n-y1nξ+y1n,0⩽ξ⩽1,where (x1n,y1n) and (x2n,y2n) are the coordinates of the ending points of the line, whose length at the end of the last time step can be computed by(39)sn-1=x2n-1-x1n-12+y2n-1-y1n-12,where time step n=1,2,…,N. Then, we can obtain the potential energy contributed by the loading {F}:(40)Πs=-∫01‍unxξ,yξvnxs,ysTFxξFyξsn-1dξ.Substituting (3) and (1) into (40) leads to(41)Πs=-DinT∫01‍TiξTFx(ξ)Fy(ξ)sn-1dξ.Minimize Πs with respect to {Dξn}, and we have(42)∂Πs0∂DinT=∫10‍TiξTFx(ξ)Fy(ξ)sn-1dξ⟶Fi.The 6×1 submatrix is added to the right item of the global equation (5). If {F}=(Fx,Fy), which is constant, (42) has the analytical form:(43)∂Πs0∂DinT=sn-1FxFy2y0-y~2Fx+x~-2x02Fyx~-2x02Fxy~-2y02Fy2y0-y~4Fx+x~-2x04Fy⟶Fi,where x~=x2+x1,y~=y2+y1 and (x0,y0) is the coordinate of the barycenter of the element Ii.


### 3.6. Volume Force Submatrix

When the elementIi bears the body loading f(x,y)=(fx(x,y), fy(x,y))T,i=1,2,…,M, the potential energy Πw due to the force f is calculated as(44)Πw=-∬Ii‍unx,yvnx,yTfx,ydxdy.Substituting (1) into the equations above leads to(45)Πw=-DinT∬Ii‍Tix,yTfx,ydxdy.Minimize Πw with respect to {Dξn}, and we have(46)∂Πw(0)∂DinT=-∬Ii‍Tix,yTfx,ydxdy⟶Fi.The 6×1 submatrix above is added to the right side of the global simultaneous equation (5). When the body force is constant, say f(x,y)=(fx,fy)T, (46) becomes(47)∂Πw(0)∂DinT=-DinTfxAi,fyAi,0,0,0,0T⟶Fi,where Ai denotes the area of the element Ii.


### 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 Figure6(a), the vertex P1 of the element Ii and the edge P2P3¯ of the element Ij 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 P0 on the edge P2P3¯ is the contact point. P1 is a vertex of the block which is different from the block where P0, P2, and P3 locate. Denote the coordinate of Pα at tn by (xαn,yαn), α=0,1,2,3 and n=1,2,…,N. The displacement of Pα in the nth time step can be expressed as(48)U(xαn,yαn)=uαnvαn=xαn-xαn-1yαn-yαn-1.Substitute (1) into the above equation, and we have(49)u1nv1n=Ti1Din,(50)uαnvαn=TjβDjn,β=0,2,3,where [Tξα] is the abbreviation of [Tξ(xα,yα)], α=0,1,2,3.Figure 6
Contact springs between the elementsIi and Ij. (a) The normal contact spring with the stiffness κ⊥ and the shear contact spring with the stiffness κ∥. (b) The normal contact distance h⊥ and the shear contact distance h∥.
(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 as(51a)Π⊥=12κ⊥h⊥2,(51b)Π∥=12κ∥h∥2,where h⊥ and h∥ denote the normal and the shear contact displacement, respectively. As shown in Figure 6(b), h⊥ is the distance from P1 to the edge P2P3¯ and can be computed as follows:(52)h⊥=AnLn,where An is the area of the triangle ΔP1P2P3 and Ln denotes the length of the edge P2P3¯. Since h∥ is the projection of the line P0P1¯ on the edge P2P3¯ or its extension line, we have(53)h∥=P2P3→·P0P1→P2P3¯=1Lnx3n-x2ny3n-y2nTx1n-x0ny1n-y0n.As 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 h⊥ firstly. An in (52) can be deduced by using (48) as(54)An=1x1ny1n1x2ny2n1x3ny3n=1x1n-1+u1ny1n-1+v1n1x2n-1+u2ny2n-1+v2n1x3n-1+u3ny3n-1+v3n=1x1n-1y1n-11x2n-1y2n-11x3n-1y3n-1+1u1ny1n-11u2ny2n-11u3ny3n-1+1x1n-1v1n1x2n-1v2n1x3n-1v3n+1u1nv1n1u2nv2n1u3nv3n.For a small enough time step, the last item in (54) will be a second-order infinitesimal and hence can be ignored. Then, we have(55)An≈An-1+1u1ny1n-11u2ny2n-11u3ny3n-1+1x1n-1v1n1x2n-1v2n1x3n-1v3n.In the same case, Ln can be approximated as Ln-1:(56)Ln≈Ln-1=x2n-1-x3n-12+y2n-1-y3n-12.Substituting (50), (55), and (56) into (52) leads to(57)h⊥=An-1Ln-1+1Ln-1y2n-1-y3n-1x3n-1-x2n-1TTi1Din+1Ln-1y3n-1-y1n-1x1n-1-x3n-1TTj2Djn+1Ln-1y1n-1-y2n-1x2n-1-x1n-1TTj3Djn.Substitute (57) into (51a), and then minimize Π⊥ with respect to {Din} and {Djn}:(58a)∂2Π⊥∂DinT∂Din=κ⊥Φ⊥TΦ⊥⟶Kii,(58b)∂2Π⊥∂DinT∂Djn=κ⊥Φ⊥TΨ⊥⟶Kij,(58c)∂2Π⊥∂DjnT∂Din=κ⊥Ψ⊥TΦ⊥⟶Kji,(58d)∂2Π⊥∂DjnT∂Djn=κ⊥Ψ⊥TΨ⊥⟶Kjj,(58e)-∂Π⊥0∂Din=κ⊥An-1Φ⊥⟶Fi,(58f)-∂Π⊥0∂Djn=κ⊥An-1Φ⊥⟶Fj,where(59a)Φ⊥=1Ln-1y2n-1-y3n-1x3n-1-x2n-1TTi1,(59b)Ψ⊥=1Ln-1y3n-1-y1n-1x1n-1-x3n-1TTj2+1Ln-1y1n-1-y2n-1x2n-1-x1n-1TTj2.The 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 h∥ first. To this end, substitute (48) into (53), and we have(60)h∥=1Lnx3n-1-x2n-1y3n-1-y2n-1Tx1n-1-x0n-1y1n-1-y0n-1+1Lnx3n-1-x2n-1y3n-1-y2n-1Tu1n-u0nv1n-v0n+1Lnu3n-u2nv3n-v2nTx1n-1-x0n-1y1n-1-y0n-1+1Lnu3n-u2nv3n-v2nTu1n-u0nv1n-v0n.With a small enough stepsize, the last two items in the equation above are infinitesimals and hence can be ignored, while Ln can be approximately replaced by Ln-1. Therefore, substituting (50) into (60) leads to(61)h∥=Bn-1Ln-1+1Ln-1x3n-1-x2n-1y3n-1-y2n-1TTi1Din+1Ln-1x2n-1-x3n-1y2n-1-y3n-1TTj0Djn,where(62)Bn-1=x1n-1-x0n-1y1n-1-y0n-1Tx3n-1-x2n-1y3n-1-y2n-1.Now, substitute (61) into (51b) and then minimize Π∥ with respect to {Din} and {Djn}, and we obtain(63a)∂2Π∥∂DinT∂Din=κ∥Φ∥⊤Φ∥⟶Kii,(63b)∂2Π∥∂DinT∂Djn=κ∥Φ∥⊤Ψ∥⟶Kij,(63c)∂2Π∥∂DjnT∂Din=κ∥Ψ∥⊤Φ∥⟶Kji,(63d)∂2Π∥∂DjnT∂Djn=κ∥Ψ∥⊤Ψ∥⟶Kjj,(63e)-∂Π∥0∂Din=κ∥Bn-1Φ∥⟶Fi,(63f)-∂Π∥0∂Djn=κ∥Bn-1Φ∥⟶Fj,which are added to the global stiffness matrix and the right item in (5), respectively, and here,(64)[Φ∥]=1Ln-1y2n-1-y3n-1,x3n-1-x2n-1Ti1,[Ψ∥]=1Ln-1x2n-1-x3n-1,y2n-1-y3n-1Tj0.


### 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 Figure7(a), the vertex P1 of the element Ii is sliding along the edge P2P3¯ of the element Ij, i,j=1,2,…,M. A normal contact spring with the stiffness κ⊥ is implemented to attempt to push the element Ii out of the element Ij. As shown in Figure 7(b), the point P0 on P2P3¯ is the contact point, s0 and s1 are the sliding displacements of P0 and P1 along the directed edge P2P3→, respectively, and h⊥ 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:(65a)N=κ⊥h⊥,(65b)F=sgnP0P1→·P2P3→Ntan⁡φ,where the sign function sgn(x) takes the values of 1, 0, and −1 in the cases of x>0, x=0, and x<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:(66)ΠF=Fs0+Fs1.At the time node tn, n=1,2,…,N, assume the coordinate of the point Pα is (xαn,yαn), α=0,1,2,3. In this way, the displacements s0 and s1 during the nth time step can be computed by(67a)s0=1Lnu0n,v0nP2P3→=1Lnu0n,v0nx3n-1-x2n-1y3n-1-y2n-1,(67b)s1=1Lnu1n,v1nP2P3→=1Lnu1n,v1nx3n-1-x2n-1y3n-1-y2n-1,where Ln denotes the length of P2P3¯ at tn. For a small enough time interval [tn-1,tn], Ln can be approximated by Ln-1, which can be computed by (56). In this case, substituting (1) and (3) into (67a) and (67b) leads to(68a)s0=1Ln-1DjnTTj0Tx3n-1-x2n-1y3n-1-y2n-1,(68b)s1=1Ln-1DinTTi1Tx3n-1-x2n-1y3n-1-y2n-1.Substitute (68a) and (68b) into (66) and then minimize ΠF with respect to {Din} and {Djn}, respectively, and we have(69a)-∂ΠF0∂Din=1Ln-1FTi1Tx2n-1-x3n-1y2n-1-y3n-1⟶Fi,(69b)-∂ΠF0∂Djn=1Ln-1FTj0Tx3n-1-x2n-1y3n-1-y2n-1⟶Fj.The above 6×1 submatrices are added to the right item of (5).Figure 7
The vertexP1 of the element Ii is sliding along the edge P2P3¯ of the element Ij. (a) A normal contact spring with the stiffness κ⊥ is applied (b).
(a)
(b)


## 3.1. Element Elastic Submatrices

The strain energy of theith element Ii can be expressed as(6)Πei=12∬Ii‍εinTσin+2σ¯indxdy,where the vectors {εin}, {σin}, and {σ¯in} denote the strain, stress, and the initial stress in the nth time step, respectively. For isotropic linearly elastic materials, {σin} and {εin} have the following relationship:(7)σin=E1-ν21ν0ν10001-ν2εin=Eiεin,in which E and ν indicate Young’s modulus and Poisson’s ratio, respectively. The strain {εin} can be expressed by the displacement as(8)εin=000100000010000001Din=ΘDin.Substituting (7) and (8) into (6), we have(9)Πei=12DinT∬Ii‍ΘTEiΘDin+DinT∬Ii‍σ¯indxdy.Minimizing Πei with respect to {Din} leads to(10)∂2Πei∂DinT∂Din=∬Ii‍ΘTEiΘdxdy⟶Kii,-∂Πei0∂DinT=-∬Ii‍ΘTσ¯indxdy⟶Fi,which 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 elementIi at the moment t can be expressed as follows:(11)Πρi=∬Ii‍Uit,x,yTρ∂2∂t2Uit,x,ydxdy,where {Ui(t,x,y)} 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 Ii. Substituting (1) into (11) leads to(12)Πρi=DitT∬Ii‍ρTiTTi∂2Dit∂t2dxdy.Using the Newmark time integrating method, we have the following finite difference approximation:(13)∂2Dit∂t2=2DinΔn2-2V¯in-1Δn,where Δn=tn-tn-1 is the nth stepsize and V¯in-1 indicates the initial velocity of the nth time step and can be updated by the iteration formulation below:(14)V¯in=2Din-1Δn-V¯in-1,n⩾2.Substituting (13) into (12), we then have(15)Πρi=2DinTΔn2∬Ii‍ρTiTTidxdyDin-2DinTΔn∬Ii‍ρTiTTidxdyV¯in-1.Minimizing the right items of the equation above with respect to {Din} leads to(16)2Δn2∬Ii‍ρTiTTidxdy⟶Kii,(17)2Δn∬Ii‍ρTiTTidxdyV¯in-1⟶Fi,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.


## 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 Figure4(a), the points P1 and P2 are a pair of initially overlapped vertexes of the elements Ii and Ij which are within the same block. A bonding spring with the stiffness κ is applied between P1 and P2 to glue Ii and Ij together. Denote the coordinate of P1 and P2 at the time node tn by (x1n,y1n) and (x2n,y2n), respectively.Figure 4
Bonding springs applied along overlapped element edges. (a) Bonding spring length is zero att0. (b) The bonding spring between P1 and P2 has the elongation dn at tn.
(a)
(b)Firstly, we compute the elongationdn of the bonding spring between P1 and P2 at the time node tn. As illustrated in Figure 4(b), dn is actually the length of P1P2→ and can be computed from the following equality:(18)dn2=lnTln,where {ln} denotes the vector P1P2→ at tn, that is,(19)ln=x2ny2n-x1ny1n.Assume the displacements of P1 and P2 in the nth time step are (u1n,v1n) and (u2n,v2n), respectively. Then, the coordinates of P1 and P2 at tn can be decomposed into the incremental form below:(20a)x1n,y1n=x1n-1+u1n,y1n-1+v1n,(20b)x2n,y2n=x2n-1+u2n,y2n-1+v2n.Substituting (20a) and (20b) into (19) leads to(21)ln=ln-1+u1nv1n-u2nv2n.Substitute (3) and (1) into the equation above, and then we have(22){ln}=ln-1+Ti1Din-Tj2Djn,where [Ti1]=[Ti(x1,y1)] and [Tj2]=[Tj(x2,y2)]. Substituting (22) into (18), then dn2 can be expanded as the following summation:(23)dn2=ln-1Tln-1+DinTTi1TTi1Din+DjnTTj2TTj2Djn-DinTTi1TTj2Djn-DjnTTj2TTi1Din+ln-1TTi1Din+Ti1TDinTln-1-ln-1TTj2Djn-Tj2TDjnTln-1.Next, we compute the strain energyΠgκ of the bonding spring between P1 and P2 in the nth time step [tn-1,tn]. Πgκ can be expressed as follows:(24)Πgκ=12κdn2.Substitute (23) into the expression above, and Πgκ becomes the following sum:(25)Πgκ=π0+πi+πj+πii+πij+πji+πjj,where(26a)π0=κ2ln-1Tln-1,(26b)πii=κ2DinTTi1TTi1Din,(26c)πjj=κ2DjnTTj2TTj2Djn,(26d)πij=-κ2DinTTi1TTj2Djn,(26e)πji=-κ2DjnTTj2TTi1Din,(26f)πi=κ2ln-1TTi1Din+κ2Ti1TDinTln-1,(26g)πj=-κ2ln-1TTj2Djn-κ2Tj2TDjnTln-1.Then, minimizeΠgκ with respect to the unknowns. As {ln-1} is known in the nth time step, the derivative of π0 is 0. Minimize πii, πij, πji, and πjj with respect to {Di} and {Dj}, and we have(27a)∂2πii∂DinT∂Din=κTi1TTi1⟶Kii,(27b)∂2πij∂DinT∂Djn=κTi1TTj2⟶Kij,(27c)∂2πji∂DjnT∂Din=κTj2TTi1⟶Kji,(27d)∂2πjj∂DjnT∂Djn=κTj2TTj2⟶Kjj.These 6×6 submatrices in the equation above are added into the global stiffness matrix in (5). Minimizing πi and πj with respect to {Di} and {Dj} leads to(28a)-∂πi0∂Di=-κTi1Tln-1⟶Fi,(28b)-∂πj({0})∂Dj=κTj2Tln-1⟶Fj,which 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 Figure5(a), the point P of the element Ii is a fixed point, i=1,2,…,M. Denote the coordinates of P at the time node tn by (xPn,yPn), n=1,2,…,N. Assume there exists a fictitious element I0 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 P0 of I0 overlaps with the fixed point P at the initial time t0. In order to resist the displacement of the point P, a bonding spring with the stiffness κ0 is applied between P and P0. The stiffness κ0 can be determined by the following formulation:(29)κ0=κ~E,where κ~ is the same factor used in (4).Figure 5
Bonding springs are applied at fixed points. (a) The length of the bonding spring applied between the fixed pointP of the element Ii and the vertex P0 of the fictitious support element I0 is zero at t0. (b) At tn, the bonding spring between P0 and P has the elongation dn.
(a)
(b)As shown in Figure5(b), the bonding spring between P and P0 has the elongation dn at the time node tn. Denote the coordinate of P0 by (x¯,y¯), which keeps constant during a computation. Since dn is actually the length of P0P→, we have(30)dn2=lPnTlPn,where {lPn} denotes the vector P0P→ at tn; that is,(31)lPn=xPnyPn-x¯y¯.Assume the displacement of P in the nth time step is (uPn,vPn). Then, the coordinates of P at tn can be written as xPn,yPn=xPn-1+uPn,yPn-1+vPn, which is substituted into (31); that is,(32)lPn=xPn-1+uPnyPn-1+vPn-x¯y¯=ln-1+uPnvPn.Substitute (3) and (1) into the equation above, and then we have(33)lPn=lPn-1+TiPDin,where [TiP] indicates [Ti(xPn,yPn)]. Substituting (33) into (30) leads to(34)dn2=lPn-1TlPn-1+DinTTi1TTi1Din+lPn-1TTi1Din+Ti1TDinTlPn-1.Then we can express the strain energy of the bonding spring between P and P0 as(35)ΠPκ0=12κ0dn2=π0+πii+πi,where(36a)π0=κ02lPn-1TlPn-1,(36b)πii=κ02DinTTi1TTi1Din,(36c)πi=κ0Ti1TDinTlPn-1.Obviously, the minimization of π0 with respect to {Din} is zero. Then, we minimize πii and πi:(37a)∂2πii∂DinT∂Din=κ0TPTTP⟶Kii,(37b)∂πi∂DinT=κ0lPn-1TTP⟶Fi,where [Kii] and {Fi} are added to the stiffness matrix and the right item of (5).


## 3.5. Line Loading Submatrices

Assume a loading{F}=(Fx(x,y),Fy(x,y))T is applied along a line upon the boundary of the element Ii, i=1,2,…,M. The coordinate of any point on the line at the time moment tn can be determined by the following parametric equations:(38)xn=xnξ=x2n-x1nξ+x1n,yn=ynξ=y2n-y1nξ+y1n,0⩽ξ⩽1,where (x1n,y1n) and (x2n,y2n) are the coordinates of the ending points of the line, whose length at the end of the last time step can be computed by(39)sn-1=x2n-1-x1n-12+y2n-1-y1n-12,where time step n=1,2,…,N. Then, we can obtain the potential energy contributed by the loading {F}:(40)Πs=-∫01‍unxξ,yξvnxs,ysTFxξFyξsn-1dξ.Substituting (3) and (1) into (40) leads to(41)Πs=-DinT∫01‍TiξTFx(ξ)Fy(ξ)sn-1dξ.Minimize Πs with respect to {Dξn}, and we have(42)∂Πs0∂DinT=∫10‍TiξTFx(ξ)Fy(ξ)sn-1dξ⟶Fi.The 6×1 submatrix is added to the right item of the global equation (5). If {F}=(Fx,Fy), which is constant, (42) has the analytical form:(43)∂Πs0∂DinT=sn-1FxFy2y0-y~2Fx+x~-2x02Fyx~-2x02Fxy~-2y02Fy2y0-y~4Fx+x~-2x04Fy⟶Fi,where x~=x2+x1,y~=y2+y1 and (x0,y0) is the coordinate of the barycenter of the element Ii.


## 3.6. Volume Force Submatrix

When the elementIi bears the body loading f(x,y)=(fx(x,y), fy(x,y))T,i=1,2,…,M, the potential energy Πw due to the force f is calculated as(44)Πw=-∬Ii‍unx,yvnx,yTfx,ydxdy.Substituting (1) into the equations above leads to(45)Πw=-DinT∬Ii‍Tix,yTfx,ydxdy.Minimize Πw with respect to {Dξn}, and we have(46)∂Πw(0)∂DinT=-∬Ii‍Tix,yTfx,ydxdy⟶Fi.The 6×1 submatrix above is added to the right side of the global simultaneous equation (5). When the body force is constant, say f(x,y)=(fx,fy)T, (46) becomes(47)∂Πw(0)∂DinT=-DinTfxAi,fyAi,0,0,0,0T⟶Fi,where Ai denotes the area of the element Ii.


## 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 Figure6(a), the vertex P1 of the element Ii and the edge P2P3¯ of the element Ij 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 P0 on the edge P2P3¯ is the contact point. P1 is a vertex of the block which is different from the block where P0, P2, and P3 locate. Denote the coordinate of Pα at tn by (xαn,yαn), α=0,1,2,3 and n=1,2,…,N. The displacement of Pα in the nth time step can be expressed as(48)U(xαn,yαn)=uαnvαn=xαn-xαn-1yαn-yαn-1.Substitute (1) into the above equation, and we have(49)u1nv1n=Ti1Din,(50)uαnvαn=TjβDjn,β=0,2,3,where [Tξα] is the abbreviation of [Tξ(xα,yα)], α=0,1,2,3.Figure 6
Contact springs between the elementsIi and Ij. (a) The normal contact spring with the stiffness κ⊥ and the shear contact spring with the stiffness κ∥. (b) The normal contact distance h⊥ and the shear contact distance h∥.
(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 as(51a)Π⊥=12κ⊥h⊥2,(51b)Π∥=12κ∥h∥2,where h⊥ and h∥ denote the normal and the shear contact displacement, respectively. As shown in Figure 6(b), h⊥ is the distance from P1 to the edge P2P3¯ and can be computed as follows:(52)h⊥=AnLn,where An is the area of the triangle ΔP1P2P3 and Ln denotes the length of the edge P2P3¯. Since h∥ is the projection of the line P0P1¯ on the edge P2P3¯ or its extension line, we have(53)h∥=P2P3→·P0P1→P2P3¯=1Lnx3n-x2ny3n-y2nTx1n-x0ny1n-y0n.As 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 h⊥ firstly. An in (52) can be deduced by using (48) as(54)An=1x1ny1n1x2ny2n1x3ny3n=1x1n-1+u1ny1n-1+v1n1x2n-1+u2ny2n-1+v2n1x3n-1+u3ny3n-1+v3n=1x1n-1y1n-11x2n-1y2n-11x3n-1y3n-1+1u1ny1n-11u2ny2n-11u3ny3n-1+1x1n-1v1n1x2n-1v2n1x3n-1v3n+1u1nv1n1u2nv2n1u3nv3n.For a small enough time step, the last item in (54) will be a second-order infinitesimal and hence can be ignored. Then, we have(55)An≈An-1+1u1ny1n-11u2ny2n-11u3ny3n-1+1x1n-1v1n1x2n-1v2n1x3n-1v3n.In the same case, Ln can be approximated as Ln-1:(56)Ln≈Ln-1=x2n-1-x3n-12+y2n-1-y3n-12.Substituting (50), (55), and (56) into (52) leads to(57)h⊥=An-1Ln-1+1Ln-1y2n-1-y3n-1x3n-1-x2n-1TTi1Din+1Ln-1y3n-1-y1n-1x1n-1-x3n-1TTj2Djn+1Ln-1y1n-1-y2n-1x2n-1-x1n-1TTj3Djn.Substitute (57) into (51a), and then minimize Π⊥ with respect to {Din} and {Djn}:(58a)∂2Π⊥∂DinT∂Din=κ⊥Φ⊥TΦ⊥⟶Kii,(58b)∂2Π⊥∂DinT∂Djn=κ⊥Φ⊥TΨ⊥⟶Kij,(58c)∂2Π⊥∂DjnT∂Din=κ⊥Ψ⊥TΦ⊥⟶Kji,(58d)∂2Π⊥∂DjnT∂Djn=κ⊥Ψ⊥TΨ⊥⟶Kjj,(58e)-∂Π⊥0∂Din=κ⊥An-1Φ⊥⟶Fi,(58f)-∂Π⊥0∂Djn=κ⊥An-1Φ⊥⟶Fj,where(59a)Φ⊥=1Ln-1y2n-1-y3n-1x3n-1-x2n-1TTi1,(59b)Ψ⊥=1Ln-1y3n-1-y1n-1x1n-1-x3n-1TTj2+1Ln-1y1n-1-y2n-1x2n-1-x1n-1TTj2.The 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 h∥ first. To this end, substitute (48) into (53), and we have(60)h∥=1Lnx3n-1-x2n-1y3n-1-y2n-1Tx1n-1-x0n-1y1n-1-y0n-1+1Lnx3n-1-x2n-1y3n-1-y2n-1Tu1n-u0nv1n-v0n+1Lnu3n-u2nv3n-v2nTx1n-1-x0n-1y1n-1-y0n-1+1Lnu3n-u2nv3n-v2nTu1n-u0nv1n-v0n.With a small enough stepsize, the last two items in the equation above are infinitesimals and hence can be ignored, while Ln can be approximately replaced by Ln-1. Therefore, substituting (50) into (60) leads to(61)h∥=Bn-1Ln-1+1Ln-1x3n-1-x2n-1y3n-1-y2n-1TTi1Din+1Ln-1x2n-1-x3n-1y2n-1-y3n-1TTj0Djn,where(62)Bn-1=x1n-1-x0n-1y1n-1-y0n-1Tx3n-1-x2n-1y3n-1-y2n-1.Now, substitute (61) into (51b) and then minimize Π∥ with respect to {Din} and {Djn}, and we obtain(63a)∂2Π∥∂DinT∂Din=κ∥Φ∥⊤Φ∥⟶Kii,(63b)∂2Π∥∂DinT∂Djn=κ∥Φ∥⊤Ψ∥⟶Kij,(63c)∂2Π∥∂DjnT∂Din=κ∥Ψ∥⊤Φ∥⟶Kji,(63d)∂2Π∥∂DjnT∂Djn=κ∥Ψ∥⊤Ψ∥⟶Kjj,(63e)-∂Π∥0∂Din=κ∥Bn-1Φ∥⟶Fi,(63f)-∂Π∥0∂Djn=κ∥Bn-1Φ∥⟶Fj,which are added to the global stiffness matrix and the right item in (5), respectively, and here,(64)[Φ∥]=1Ln-1y2n-1-y3n-1,x3n-1-x2n-1Ti1,[Ψ∥]=1Ln-1x2n-1-x3n-1,y2n-1-y3n-1Tj0.


## 3.7.1. Normal Contact Submatrices

To obtain the derivative ofΠ⊥ in (51a), we compute the normal contact displacement h⊥ firstly. An in (52) can be deduced by using (48) as(54)An=1x1ny1n1x2ny2n1x3ny3n=1x1n-1+u1ny1n-1+v1n1x2n-1+u2ny2n-1+v2n1x3n-1+u3ny3n-1+v3n=1x1n-1y1n-11x2n-1y2n-11x3n-1y3n-1+1u1ny1n-11u2ny2n-11u3ny3n-1+1x1n-1v1n1x2n-1v2n1x3n-1v3n+1u1nv1n1u2nv2n1u3nv3n.For a small enough time step, the last item in (54) will be a second-order infinitesimal and hence can be ignored. Then, we have(55)An≈An-1+1u1ny1n-11u2ny2n-11u3ny3n-1+1x1n-1v1n1x2n-1v2n1x3n-1v3n.In the same case, Ln can be approximated as Ln-1:(56)Ln≈Ln-1=x2n-1-x3n-12+y2n-1-y3n-12.Substituting (50), (55), and (56) into (52) leads to(57)h⊥=An-1Ln-1+1Ln-1y2n-1-y3n-1x3n-1-x2n-1TTi1Din+1Ln-1y3n-1-y1n-1x1n-1-x3n-1TTj2Djn+1Ln-1y1n-1-y2n-1x2n-1-x1n-1TTj3Djn.Substitute (57) into (51a), and then minimize Π⊥ with respect to {Din} and {Djn}:(58a)∂2Π⊥∂DinT∂Din=κ⊥Φ⊥TΦ⊥⟶Kii,(58b)∂2Π⊥∂DinT∂Djn=κ⊥Φ⊥TΨ⊥⟶Kij,(58c)∂2Π⊥∂DjnT∂Din=κ⊥Ψ⊥TΦ⊥⟶Kji,(58d)∂2Π⊥∂DjnT∂Djn=κ⊥Ψ⊥TΨ⊥⟶Kjj,(58e)-∂Π⊥0∂Din=κ⊥An-1Φ⊥⟶Fi,(58f)-∂Π⊥0∂Djn=κ⊥An-1Φ⊥⟶Fj,where(59a)Φ⊥=1Ln-1y2n-1-y3n-1x3n-1-x2n-1TTi1,(59b)Ψ⊥=1Ln-1y3n-1-y1n-1x1n-1-x3n-1TTj2+1Ln-1y1n-1-y2n-1x2n-1-x1n-1TTj2.The 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 h∥ first. To this end, substitute (48) into (53), and we have(60)h∥=1Lnx3n-1-x2n-1y3n-1-y2n-1Tx1n-1-x0n-1y1n-1-y0n-1+1Lnx3n-1-x2n-1y3n-1-y2n-1Tu1n-u0nv1n-v0n+1Lnu3n-u2nv3n-v2nTx1n-1-x0n-1y1n-1-y0n-1+1Lnu3n-u2nv3n-v2nTu1n-u0nv1n-v0n.With a small enough stepsize, the last two items in the equation above are infinitesimals and hence can be ignored, while Ln can be approximately replaced by Ln-1. Therefore, substituting (50) into (60) leads to(61)h∥=Bn-1Ln-1+1Ln-1x3n-1-x2n-1y3n-1-y2n-1TTi1Din+1Ln-1x2n-1-x3n-1y2n-1-y3n-1TTj0Djn,where(62)Bn-1=x1n-1-x0n-1y1n-1-y0n-1Tx3n-1-x2n-1y3n-1-y2n-1.Now, substitute (61) into (51b) and then minimize Π∥ with respect to {Din} and {Djn}, and we obtain(63a)∂2Π∥∂DinT∂Din=κ∥Φ∥⊤Φ∥⟶Kii,(63b)∂2Π∥∂DinT∂Djn=κ∥Φ∥⊤Ψ∥⟶Kij,(63c)∂2Π∥∂DjnT∂Din=κ∥Ψ∥⊤Φ∥⟶Kji,(63d)∂2Π∥∂DjnT∂Djn=κ∥Ψ∥⊤Ψ∥⟶Kjj,(63e)-∂Π∥0∂Din=κ∥Bn-1Φ∥⟶Fi,(63f)-∂Π∥0∂Djn=κ∥Bn-1Φ∥⟶Fj,which are added to the global stiffness matrix and the right item in (5), respectively, and here,(64)[Φ∥]=1Ln-1y2n-1-y3n-1,x3n-1-x2n-1Ti1,[Ψ∥]=1Ln-1x2n-1-x3n-1,y2n-1-y3n-1Tj0.


## 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 Figure7(a), the vertex P1 of the element Ii is sliding along the edge P2P3¯ of the element Ij, i,j=1,2,…,M. A normal contact spring with the stiffness κ⊥ is implemented to attempt to push the element Ii out of the element Ij. As shown in Figure 7(b), the point P0 on P2P3¯ is the contact point, s0 and s1 are the sliding displacements of P0 and P1 along the directed edge P2P3→, respectively, and h⊥ 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:(65a)N=κ⊥h⊥,(65b)F=sgnP0P1→·P2P3→Ntan⁡φ,where the sign function sgn(x) takes the values of 1, 0, and −1 in the cases of x>0, x=0, and x<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:(66)ΠF=Fs0+Fs1.At the time node tn, n=1,2,…,N, assume the coordinate of the point Pα is (xαn,yαn), α=0,1,2,3. In this way, the displacements s0 and s1 during the nth time step can be computed by(67a)s0=1Lnu0n,v0nP2P3→=1Lnu0n,v0nx3n-1-x2n-1y3n-1-y2n-1,(67b)s1=1Lnu1n,v1nP2P3→=1Lnu1n,v1nx3n-1-x2n-1y3n-1-y2n-1,where Ln denotes the length of P2P3¯ at tn. For a small enough time interval [tn-1,tn], Ln can be approximated by Ln-1, which can be computed by (56). In this case, substituting (1) and (3) into (67a) and (67b) leads to(68a)s0=1Ln-1DjnTTj0Tx3n-1-x2n-1y3n-1-y2n-1,(68b)s1=1Ln-1DinTTi1Tx3n-1-x2n-1y3n-1-y2n-1.Substitute (68a) and (68b) into (66) and then minimize ΠF with respect to {Din} and {Djn}, respectively, and we have(69a)-∂ΠF0∂Din=1Ln-1FTi1Tx2n-1-x3n-1y2n-1-y3n-1⟶Fi,(69b)-∂ΠF0∂Djn=1Ln-1FTj0Tx3n-1-x2n-1y3n-1-y2n-1⟶Fj.The above 6×1 submatrices are added to the right item of (5).Figure 7
The vertexP1 of the element Ii is sliding along the edge P2P3¯ of the element Ij. (a) A normal contact spring with the stiffness κ⊥ is applied (b).
(a)
(b)


## 4. Numerical Examples

The computer program for the BBM has been coded in the C language, whose flowchart is illustrated in Figure8. 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.Figure 8
The flowchart of the computer program for the BBM.


### 4.1. Parameter Settings

As shown in Table1, the parameters used in the experiments are given, including time step interval Δt, total time steps N, the error tolerance ϵt 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 κ∥.Table 1
The values of the parameters used in the experiments.
Number
Δ
t
N
ϵ
t
ω
ϵ
κ
~
κ
⊥
κ
∥
I
1.0
100
1.0
1.0
1
e
-
8
10
—
—
II
1.0
1
1.0
1.0
1
e
-
8
1
0
2
—
—
III
1
e
-
6
40
1.0
1.0
1
e
-
8
10
1
0
2
E
†
1
0
2
E
†
E 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 Figure9(a), a cantilever beam subjected to a concentrated load P=1 KN at the free end is considered. The length and the height of the beam are L=10m and h=1m, respectively. The material parameters of the beam are as follows: Young’s modulus E=300MPa, Poisson’s ratio ν=0.3, and the mass density ρ=0.01Kg/m2. The displacement of this problem has the following analytical solution [36]: (70)Ux=Py6EI6L-3xx+2+νy2-h24,Uy=Py6EI3νy2(L-x)+(4+5ν)h2y4+(3L-x)x2,where I=h3/12 is the moment of inertia.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.
(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 Figure9(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.Table 2
Comparing the DOFs involved in the DDA, the FEM, and the BBM for Section4.2.
Models
DDA
FEM
BBM
Entity
Blocks
Element
BBM element
DOFs/entity
6
12
6
Total entities
1
2202
2202
Total DOFs
6
26424
13202This 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 numberN or to satisfy the following inequality constraint:(71)∑(x,y)∈∪i=1MIi‍Un-Un-1221/2⩽ϵNR,where time step n=1,2,…,100; element number M=2202; {Un} = {Un(x,y)} denotes the computed displacement at (x,y) in the nth time step; ·2 is the square norm; ϵNR=10-8 is the error tolerance for Newton-Raphson iterations. In this example, the computation of the BBM converges at n=3.Comparative Task 1. Figure 10 plots the computed displacements on the nodes distributed along the line AB¯ = {(x,y):0⩽x⩽L,y=-h/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.Figure 10
The computed displacement componentsUx and Uy of the nodes along the line AB.
(a)
(b)Comparative Task 2. In Table 3, the displacement components Ux and Uy at the point B(L,-h/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.Table 3
The displacements atB(L,-h/2) obtained by the analytical solution, the BBM, and the FEM.
Displ./m
Theory
BBM
Accuracy
FEM
Accuracy
U
x
(
B
)
−0.001
−0.000993
99.30%
−0.000998
99.80%
U
y
(
B
)
−0.013425
−0.013417
99.94%
−0.013410
99.89%Comparative Task 3. Figures 11 and 12 plot the computed contours of the displacement components Ux and Uy 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.Figure 11
The contour plots of the displacementUx (m) throughout the beam. (a) BBM result. (b) FEM result.
(a)
(b)Figure 12
The contour plots of the displacementUy (m) throughout the beam. (a) BBM result. (b) FEM result.
(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.Figure 13
The contour plots of the von-Mises stress (Pa) throughout the beam. (a) BBM result. (b) FEM result.
(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 Figure14(a), a thin plate with a circular hole at the center is pressured by a uniform horizontal stress σ∞=1MPa. 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/m2, respectively. The radium of the hole is 2a; the length and the width of the plate are 2L and 2W, respectively, where a=0.5m, L=10m, and W=5m. Set the origin of the Cartesian coordinate system at the point o, the center of the hole, and the ox-axis toward the right. As shown in Figure 14(a), we also use the polar coordinate system (r,θ), where the pole locates at the origin o, θ is the polar angle in radian, and r is the distance in meter away from o. In the case of infinite L and W, this problem has the following analytical solution derived in [37, 38]:(72a)σr=σ∞21-a2r2[1+1-3a2r2cos⁡⁡2θ],(72b)σθ=σ∞21+a2r2-σ∞21-3a2r2cos⁡⁡2θ.Figure 14
A finite-width plate with a central circular hole subjected to horizontal uniform tensions. (a) The geometry configuration and the mechanical constants of the plate. (b) The mesh used in the computations of the BBM and the FEM.
(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 areaABCDE, is analyzed due to the symmetry, and the symmetrical boundary constraints are applied along the AE and BC, 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 A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5). According to this table, both the BBM and the FEM can give accurate enough horizontal stresses σx at A, C, D, and E, while σx(B) 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 L and W are assumed to be infinite in the analytical formulations in (72a) and (72b), but they have the specific values of L=10 and W=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.Table 4
The stressesσx computed by the BBM and the FEM are compared with the analytical solutions at the points A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5).
Stress/MPa
Theory
BBM
Accuracy
FEM
Accuracy
σ
x
(
A
)
3.0
2.949130
98.304%
3.03450
98.850%
σ
x
(
B
)
0.0
-
2.61091
e
-
2
—
-
1.13856
e
-
3
—
σ
x
(
C
)
0.993759
0.999961
99.376%
1.000023
99.370%
σ
x
(
D
)
0.998758
0.999944
99.881%
0.999960
99.880%
σ
x
(
E
)
1.005150
0.989241
98.417%
0.985517
98.047%Comparative Task 2. In Table 5, the displacement components Ux and Uy computed by the BBM and the FEM are compared at the five measurement points: A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5). 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 Ux(A), Uy(B), Uy(C), and Ux(E) 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 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 Ux(A), Uy(B), Uy(C), and Ux(E) are acceptable.Table 5
The displacementsUx and Uy computed by the BBM and the FEM are compared at the points A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5).
Displ./m
FEM
BBM
Relative error†
U
x
(
A
)
0.0
1.69968
e
-
11
—
U
y
(
A
)
-
1.78046
e
-
07
-
1.75791
e
-
07
1.266527%
U
x
(
B
)
5.23535
e
-
07
5.22699
e
-
07
0.159684%
U
y
(
B
)
0.0
-
3.2779
e
-
12
—
U
x
(
C
)
3.49394
e
-
06
3.49639
e
-
06
0.070121%
U
y
(
C
)
0.0
9.402
e
-
13
—
U
x
(
D
)
3.48394
e
-
06
3.4873
e
-
06
0.096443%
U
y
(
D
)
-
5.09281
e
-
07
-
5.09625
e
-
07
0.067546%
U
x
(
E
)
0.0
6.93519
e
-
11
—
U
y
(
E
)
-
5.63371
e
-
07
-
5.6286624
e
-
07
0.089596%
†The relative error between the BBM result UαFEM and the FEM solution UαFEM is determined by UαFEM-UαBBM/UFEM, α=x,y.Comparative Task 3. Then, we compare the computed deformations along the ox- and oy-axis with the theoretical results. Figure 15 plots the displacements ux and uy obtained from the BBM and the FEM results on the nodes along the boundary BC¯. And, in Figure 16, ux and uy solved by the two methods on the nodes along AE¯ are drawn. From the figures, it can be observed that the BBM results are in accord with the FEM solutions.Figure 15
The displacementsUx (m) and Uy (m) computed by the BBM and the FEM on the mesh nodes along the symmetric line BC¯={(x,y):0.5⩽x⩽5.0,y=0}.Figure 16
The displacementsUx (m) and Uy (m) computed by the BBM, the FEM and the analytical formulations on the mesh nodes along the symmetric line AE¯={(x,y):x=0,0.5⩽y⩽5.0}.Comparative Task 4. As shown in Figures 17, 18, and 19, the computed horizontal stresses σx on the nodes along BC¯, AE¯, and the arc BA 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.Figure 17
The horizontal stressσx (Pa) computed by the BBM, the FEM and the analytical formulations in (72a) and (72b) on the mesh nodes along the symmetric line BC¯={(x,y):0.5⩽x⩽5.0,y=0}.Figure 18
The horizontal stressσx (Pa) computed by the BBM, the FEM and the analytical formulations in (72a) and (72b) on the mesh nodes along the symmetric line AE¯={(x,y):x=0,0.5⩽y⩽5.0}.Figure 19
The horizontal stressσx (Pa) computed by the BBM, the FEM and the analytical formulations in (72a) and (72b) on the mesh nodes along the arc BA={(r,θ):r=0.5,0⩽θ⩽π/2}.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 σx throughout the domain are illustrated, and it can be observed that the BBM results agree with the FEM very well.Figure 20
The contour plots of the computed distributions of the horizontal stressσx (Pa) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(a)
(b)Figure 21
The contour plots of the computed distributions of the horizontal normal stressσy (Pa) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(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.Figure 22
The contour plots of the computed distributions of the displacement componentUx (m) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(a)
(b)Figure 23
The contour plots of the computed distributions of the displacement componentUy (m) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(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 Figure24(a), a small rectangular block subjected to a horizontal traction force P slides along the top surface of a fixed flat plate. The origin of the coordinate system locates at the plate centroid, and the x 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).Figure 24
A rectangular block subjected to the tractionP slides along the surface of a fixed flat plate. (a) The geometry configuration and the mechanical parameters. (b) The subdivision.
(a)
(b)We implement the BBM to solve this problem by using the quadrilateral mesh as plotted in Figure24(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 A(-1,0.3) as marked in Figure 24(a) is concerned. The analytical solution (UxAnalyt,UyAnalyt) can be deduced based on kinematics theories; that is, UxAnalyt=50000·t2 and UyAnalyt=0. As shown in Figure 25, the computed displacement components Ux and Uy 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
The computed displacements measured at the pointA in the x and y directions, respectively.


## 4.1. Parameter Settings

As shown in Table1, the parameters used in the experiments are given, including time step interval Δt, total time steps N, the error tolerance ϵt 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 κ∥.Table 1
The values of the parameters used in the experiments.
Number
Δ
t
N
ϵ
t
ω
ϵ
κ
~
κ
⊥
κ
∥
I
1.0
100
1.0
1.0
1
e
-
8
10
—
—
II
1.0
1
1.0
1.0
1
e
-
8
1
0
2
—
—
III
1
e
-
6
40
1.0
1.0
1
e
-
8
10
1
0
2
E
†
1
0
2
E
†
E 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 Figure9(a), a cantilever beam subjected to a concentrated load P=1 KN at the free end is considered. The length and the height of the beam are L=10m and h=1m, respectively. The material parameters of the beam are as follows: Young’s modulus E=300MPa, Poisson’s ratio ν=0.3, and the mass density ρ=0.01Kg/m2. The displacement of this problem has the following analytical solution [36]: (70)Ux=Py6EI6L-3xx+2+νy2-h24,Uy=Py6EI3νy2(L-x)+(4+5ν)h2y4+(3L-x)x2,where I=h3/12 is the moment of inertia.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.
(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 Figure9(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.Table 2
Comparing the DOFs involved in the DDA, the FEM, and the BBM for Section4.2.
Models
DDA
FEM
BBM
Entity
Blocks
Element
BBM element
DOFs/entity
6
12
6
Total entities
1
2202
2202
Total DOFs
6
26424
13202This 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 numberN or to satisfy the following inequality constraint:(71)∑(x,y)∈∪i=1MIi‍Un-Un-1221/2⩽ϵNR,where time step n=1,2,…,100; element number M=2202; {Un} = {Un(x,y)} denotes the computed displacement at (x,y) in the nth time step; ·2 is the square norm; ϵNR=10-8 is the error tolerance for Newton-Raphson iterations. In this example, the computation of the BBM converges at n=3.Comparative Task 1. Figure 10 plots the computed displacements on the nodes distributed along the line AB¯ = {(x,y):0⩽x⩽L,y=-h/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.Figure 10
The computed displacement componentsUx and Uy of the nodes along the line AB.
(a)
(b)Comparative Task 2. In Table 3, the displacement components Ux and Uy at the point B(L,-h/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.Table 3
The displacements atB(L,-h/2) obtained by the analytical solution, the BBM, and the FEM.
Displ./m
Theory
BBM
Accuracy
FEM
Accuracy
U
x
(
B
)
−0.001
−0.000993
99.30%
−0.000998
99.80%
U
y
(
B
)
−0.013425
−0.013417
99.94%
−0.013410
99.89%Comparative Task 3. Figures 11 and 12 plot the computed contours of the displacement components Ux and Uy 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.Figure 11
The contour plots of the displacementUx (m) throughout the beam. (a) BBM result. (b) FEM result.
(a)
(b)Figure 12
The contour plots of the displacementUy (m) throughout the beam. (a) BBM result. (b) FEM result.
(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.Figure 13
The contour plots of the von-Mises stress (Pa) throughout the beam. (a) BBM result. (b) FEM result.
(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 Figure14(a), a thin plate with a circular hole at the center is pressured by a uniform horizontal stress σ∞=1MPa. 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/m2, respectively. The radium of the hole is 2a; the length and the width of the plate are 2L and 2W, respectively, where a=0.5m, L=10m, and W=5m. Set the origin of the Cartesian coordinate system at the point o, the center of the hole, and the ox-axis toward the right. As shown in Figure 14(a), we also use the polar coordinate system (r,θ), where the pole locates at the origin o, θ is the polar angle in radian, and r is the distance in meter away from o. In the case of infinite L and W, this problem has the following analytical solution derived in [37, 38]:(72a)σr=σ∞21-a2r2[1+1-3a2r2cos⁡⁡2θ],(72b)σθ=σ∞21+a2r2-σ∞21-3a2r2cos⁡⁡2θ.Figure 14
A finite-width plate with a central circular hole subjected to horizontal uniform tensions. (a) The geometry configuration and the mechanical constants of the plate. (b) The mesh used in the computations of the BBM and the FEM.
(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 areaABCDE, is analyzed due to the symmetry, and the symmetrical boundary constraints are applied along the AE and BC, 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 A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5). According to this table, both the BBM and the FEM can give accurate enough horizontal stresses σx at A, C, D, and E, while σx(B) 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 L and W are assumed to be infinite in the analytical formulations in (72a) and (72b), but they have the specific values of L=10 and W=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.Table 4
The stressesσx computed by the BBM and the FEM are compared with the analytical solutions at the points A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5).
Stress/MPa
Theory
BBM
Accuracy
FEM
Accuracy
σ
x
(
A
)
3.0
2.949130
98.304%
3.03450
98.850%
σ
x
(
B
)
0.0
-
2.61091
e
-
2
—
-
1.13856
e
-
3
—
σ
x
(
C
)
0.993759
0.999961
99.376%
1.000023
99.370%
σ
x
(
D
)
0.998758
0.999944
99.881%
0.999960
99.880%
σ
x
(
E
)
1.005150
0.989241
98.417%
0.985517
98.047%Comparative Task 2. In Table 5, the displacement components Ux and Uy computed by the BBM and the FEM are compared at the five measurement points: A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5). 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 Ux(A), Uy(B), Uy(C), and Ux(E) 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 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 Ux(A), Uy(B), Uy(C), and Ux(E) are acceptable.Table 5
The displacementsUx and Uy computed by the BBM and the FEM are compared at the points A(0,0.5), B(0.5,0), C(10,0), D(10,5), and E(0,5).
Displ./m
FEM
BBM
Relative error†
U
x
(
A
)
0.0
1.69968
e
-
11
—
U
y
(
A
)
-
1.78046
e
-
07
-
1.75791
e
-
07
1.266527%
U
x
(
B
)
5.23535
e
-
07
5.22699
e
-
07
0.159684%
U
y
(
B
)
0.0
-
3.2779
e
-
12
—
U
x
(
C
)
3.49394
e
-
06
3.49639
e
-
06
0.070121%
U
y
(
C
)
0.0
9.402
e
-
13
—
U
x
(
D
)
3.48394
e
-
06
3.4873
e
-
06
0.096443%
U
y
(
D
)
-
5.09281
e
-
07
-
5.09625
e
-
07
0.067546%
U
x
(
E
)
0.0
6.93519
e
-
11
—
U
y
(
E
)
-
5.63371
e
-
07
-
5.6286624
e
-
07
0.089596%
†The relative error between the BBM result UαFEM and the FEM solution UαFEM is determined by UαFEM-UαBBM/UFEM, α=x,y.Comparative Task 3. Then, we compare the computed deformations along the ox- and oy-axis with the theoretical results. Figure 15 plots the displacements ux and uy obtained from the BBM and the FEM results on the nodes along the boundary BC¯. And, in Figure 16, ux and uy solved by the two methods on the nodes along AE¯ are drawn. From the figures, it can be observed that the BBM results are in accord with the FEM solutions.Figure 15
The displacementsUx (m) and Uy (m) computed by the BBM and the FEM on the mesh nodes along the symmetric line BC¯={(x,y):0.5⩽x⩽5.0,y=0}.Figure 16
The displacementsUx (m) and Uy (m) computed by the BBM, the FEM and the analytical formulations on the mesh nodes along the symmetric line AE¯={(x,y):x=0,0.5⩽y⩽5.0}.Comparative Task 4. As shown in Figures 17, 18, and 19, the computed horizontal stresses σx on the nodes along BC¯, AE¯, and the arc BA 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.Figure 17
The horizontal stressσx (Pa) computed by the BBM, the FEM and the analytical formulations in (72a) and (72b) on the mesh nodes along the symmetric line BC¯={(x,y):0.5⩽x⩽5.0,y=0}.Figure 18
The horizontal stressσx (Pa) computed by the BBM, the FEM and the analytical formulations in (72a) and (72b) on the mesh nodes along the symmetric line AE¯={(x,y):x=0,0.5⩽y⩽5.0}.Figure 19
The horizontal stressσx (Pa) computed by the BBM, the FEM and the analytical formulations in (72a) and (72b) on the mesh nodes along the arc BA={(r,θ):r=0.5,0⩽θ⩽π/2}.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 σx throughout the domain are illustrated, and it can be observed that the BBM results agree with the FEM very well.Figure 20
The contour plots of the computed distributions of the horizontal stressσx (Pa) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(a)
(b)Figure 21
The contour plots of the computed distributions of the horizontal normal stressσy (Pa) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(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.Figure 22
The contour plots of the computed distributions of the displacement componentUx (m) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(a)
(b)Figure 23
The contour plots of the computed distributions of the displacement componentUy (m) throughout the plate. (a) The result of the BBM. (b) The result of the FEM.
(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 Figure24(a), a small rectangular block subjected to a horizontal traction force P slides along the top surface of a fixed flat plate. The origin of the coordinate system locates at the plate centroid, and the x 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).Figure 24
A rectangular block subjected to the tractionP slides along the surface of a fixed flat plate. (a) The geometry configuration and the mechanical parameters. (b) The subdivision.
(a)
(b)We implement the BBM to solve this problem by using the quadrilateral mesh as plotted in Figure24(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 A(-1,0.3) as marked in Figure 24(a) is concerned. The analytical solution (UxAnalyt,UyAnalyt) can be deduced based on kinematics theories; that is, UxAnalyt=50000·t2 and UyAnalyt=0. As shown in Figure 25, the computed displacement components Ux and Uy 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
The computed displacements measured at the pointA in the x and y directions, respectively.


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



---

*Source: 723263-2015-04-07.xml*