# Study on the Damping of the Discontinuous Deformation Analysis Based on Two-Block Model

Authors: Xixia Feng; Qinghui Jiang; Aiqing Wu; Wei Wei; Xiaobo Zhang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The deformation and failure of rock mass is a process of energy dissipation; the damping of DDA is a very important and basic problem. The correctness and effectiveness of DDA rely on the appropriate values of the numeric controlling parameters like time interval, spring stiffness, and assumed maximum displacement ratiog2, and the contact using the penalty method is the core content of DDA. A mechanical model of two contact blocks loaded with the normal force acting along one side of block boundary is established to study the DDA damping problem, which involves the contact and eliminates the influence of some numeric control parameters (e.g., g2). Based on the Newmark method and the theory of DDA, the motion equations of two-block system can be established, and then the relationship of some numeric control parameters and the influence of damping can be obtained. The algorithmic damping increases with the increasing of time interval. Given a very small time interval, the spring stiffness may have no obvious effect on the algorithmic damping. The numerical results reveal that the essence of time interval influencing the open-close iteration is the fact that the algorithmic damping is mainly controlled by time interval.


---

## Body

## 1. Introduction

The strong earthquakes can destroy rock mass and induce an abundance of landslides, which cause potentially serious threat to both lives and properties. So it is very necessary to study the seismic response of rock mass. Generally, the numerical methods are widely used for dealing with the problem. Furthermore, as the unfavorable geological structure of rock masses, the discontinue-based methods are more suitable for studying the problem.The discontinuous deformation analysis (DDA), proposed by Shi and Goodman [1, 2], can be used to simulate the discontinuous deformation behavior of jointed rock masses. The correctness and effectiveness of DDA have been verified by MacLaughlin and Doolin [3]. The DDA method is established on the block kinetics; the block dynamic analysis of DDA has been confirmed by many scholars [4–6]. Meanwhile, the DDA method has been widely applied to the seismic response of rock masses. Zhang et al. [7] investigated the seismic response analysis of large underground caverns with the modified DDA method. Zhang et al. [8] studied the run-out analysis of the Daguangbao landslides subjected to near-fault multidirection earthquake forces. Zhang et al. [9] have verified the mobility of earthquake-induced landslide by DDA method. Wu and Chen [10] simulated the kinematic behavior of sliding blocks of rock in the earthquake-induced Tsaoling landslide using seismic discontinuous deformation analysis (DDA). These studies show that DDA method has great potential in the analysis of seismic response of rock masses.It is generally recognized that the deformation and failure of rock mass is a process of energy dissipation, which should be highlighted in the dynamic problem of rock mass, such as earthquakes. Since the damping research of DDA is the basis of the seismic response, the damping characteristics should be investigated when using DDA to simulate the earthquake-induced dynamic problems of rock mass. In original DDA program, kinetic damping was introduced to discount the initial velocity of each time step. In static analysis the initial velocity is set to be zero for each time step, while the initial velocity for the dynamic analyses has not been discounted. It should be noted that the value of kinetic damping plays a key role in the dynamic analysis of rock mass. Hatzor et al. [11] and Tsesarsky et al. [5] obtained that 2% kinetic damping was the correct number for the dynamic analysis of a single block on an inclined subjected to dynamic loading. However, the significance of the kinetic damping is still not clear, and the value of initial velocity discounted is hard to choose; therefore, some researchers have introduced self-adaptive damping [12] and viscous damping [13, 14] to address these problems. Lin and Xie [15] also analyzed the performance of Newmark time integration with kinetic damping.Algorithmic damping is caused by the integration scheme of inertia force, which is a controlling factor of the convergence of open-close iteration. Doolin and Sitar [16] studied the time integration of DDA and noted that the algorithmic damping may be important considering the penalty formulation. Bao et al. [17] obtained an equivalent damping ratio from the algorithmic damping in DDA using the cantilever beam model. Overall, the damping researches of DDA are always in progress.It is well known that the numeric controlling parameters (e.g., assumed maximum displacement ratiog2, time interval) in DDA are hard to choose; the correctness and effectiveness of DDA depend on the appropriate values of the these numeric control parameters. Moreover, the contact using the penalty method is the core content of DDA, and the open-close iteration convergence is the highlight. Over the years, researches of DDA damping that involve the contact and eliminate the influence of some numeric control parameters (e.g., g2) have been rarely reported. In this study, a mechanical model of two contact blocks loaded with the normal force acting along one side of block boundary is established to study the aforementioned problems. Based on the mechanical model, the total equation motion is established easily and intuitively, and the influence of damping is analyzed. The results show that the essence of time interval influencing the open-close iteration is the algorithmic damping mainly controlled by the time interval.


## 2. Equations of Motion of a Discrete Block System

The DDA method is in essence a numerical method to solve the motion and deformation of multiple rock blocks. Starting from the definitions of displacement functions, the equations of motion of the blocks are established firstly using potential energy minimization. The direct time-integration scheme is adopted to solve the equations of motion.


### 2.1. Displacement Functions of Deformable Blocks

The DDA method selects the rigid body displacements and the strains of block elements as the basic unknown variables to solve for solutions. Assuming that an arbitrarily shaped polygonal block has uniform stress and strain, an incremental expression for the unknown variables of the blocki is written as(1)ΔDi=Δu0Δv0Δr0ΔεxΔεyΔγxyT,where (Δu0,Δv0) are the translation displacement increments at the point (x0,y0) in x and y directions, Δr0 is the rotation angle increment of the block around its mass center (x0,y0), and (Δεx,Δεy,Δγxy) are the increments of normal and shear strains of the block.Displacement increments(Δu,Δv) of an arbitrary point (x,y) in the block can be represented as(2)ΔuΔv=TiΔDi,where(3)Ti=10-y-y0x-x00y-y0201x-x00y-y0x-x02.


### 2.2. Equations of Motion of the Block System Based on the Newmark Method

For a discrete model, the total potential energy of the block system is the sum of all individual blocks’ potential energy. According to the minimum potential principle, the equation of motion of the block system in matrix form [13] can be expressed as(4)MD¨+KΔD=F,where D¨ is the acceleration vector and M, K, and F are the global mass matrix, stiffness matrix, and load matrix, respectively.A direct time-integration scheme is used for solving (4), in which the acceleration and the velocity within each time step are taken as the following equations of the Newmark assumption [18]:(5)D¨n+1=1αΔt2Dn+1-Dn-1αΔtD˙n-12α-1D¨n,D˙n+1=δαΔtDn+1-Dn+1-δαD˙n+1-δ2αΔtD¨n.While denoting V0=D˙(0), A0=D¨(0), and ΔD=Dn+1-Dn, it can be regarded that V0, A0, and ΔD are the initial velocity vector, the initial acceleration vector, and the displacement increment vector of the block system at the beginning of the current time step, respectively. The kinetic energy dissipation of DDA block system depends on the value of gg  (gg=0~1), which is the dynamic controlling parameter. Therefore, gg is introduced to the Newmark assumption; we have(6)D¨n+1=1αΔt2ΔD-gg1αΔtV0-12α-1A0,(7)D˙n+1=δαΔtΔD+gg1-δαV0+1-δ2αΔtA0.Substituting (6) into (4), the motion equation of the block system of the current time step can be derived as(8)K~ΔD=F~,where (9)K~=K+MαΔt2,F~=F+Mgg1αΔtV0+12α-1A0.


### 2.3. Algorithmic Damping of the Newmark Method

For analysis of the stability and accuracy of Newmark method, takingF=0, the following equation can be obtained based on (4)-(5): (10)1+αΩ2Dn+1=2-12-2α+δΩ2Dn-1+12+α-δΩ2Dn-1,where sampling frequency Ω=Δtw and undamped frequency of vibration w=K/M.The solutionsλ of the characteristic equation about (10) are(11)λ1,2=2-g±2-g2-41+h2,where (12)g=1/2+δΩ21+αΩ2,h=1/2-δΩ21+αΩ2.The real solution must have oscillation in the case of small damping, and the stability must be not infinitely growing; thus(13)41+h>2-g2,λ=1+h≤1when Δt is not restricted, (13) should be met, and then α≥(1/4)1/2+δ2, δ≥1/2, 1/2-δ+α≥0.So the conditions of the unconditionally stable Newmark method are(14)δ≥12,α≥1412+δ2.It is thus obvious that [13, 14, 16, 18] both the average acceleration method (α=1/4, δ=1/2) and the constant acceleration method (α=1/2, δ=1) meet the unconditionally stable condition, while the linear acceleration method (α=1/6, δ=1/2) and the central difference method (α=0, δ=1/2) are the conditionally stable method.In general, the accuracy of Newmark method depends on the time interval, the physical parameters of the system, and the loading condition. Here, we only analyze some widely used unconditionally stable integration schemes, to obtain an idea of the accuracy of these methods evaluated based on (11) λ=1+h. Figure 1 shows the algorithmic damping of the Newmark method versus sampling frequency Ω.Figure 1
Accuracy analysis of two Newmark integration schemes.As is shown in Figure1, λ of the average acceleration method is unity, which states that this method has no algorithmic damping. By contrast, λ of the constant acceleration method decreases with the increasing of sampling frequency Ω, which indicates that the constant acceleration method has a rapid attenuation under high frequency interference, while it has very little effect on low frequency. In fact, time interval used in Newmark method is greater than the period of the highest natural frequency of the system; thus, the frequency response is not reliable. So the algorithmic damping is advantageous under certain condition; it may be more suitable for simulating the engineering problems.


## 2.1. Displacement Functions of Deformable Blocks

The DDA method selects the rigid body displacements and the strains of block elements as the basic unknown variables to solve for solutions. Assuming that an arbitrarily shaped polygonal block has uniform stress and strain, an incremental expression for the unknown variables of the blocki is written as(1)ΔDi=Δu0Δv0Δr0ΔεxΔεyΔγxyT,where (Δu0,Δv0) are the translation displacement increments at the point (x0,y0) in x and y directions, Δr0 is the rotation angle increment of the block around its mass center (x0,y0), and (Δεx,Δεy,Δγxy) are the increments of normal and shear strains of the block.Displacement increments(Δu,Δv) of an arbitrary point (x,y) in the block can be represented as(2)ΔuΔv=TiΔDi,where(3)Ti=10-y-y0x-x00y-y0201x-x00y-y0x-x02.


## 2.2. Equations of Motion of the Block System Based on the Newmark Method

For a discrete model, the total potential energy of the block system is the sum of all individual blocks’ potential energy. According to the minimum potential principle, the equation of motion of the block system in matrix form [13] can be expressed as(4)MD¨+KΔD=F,where D¨ is the acceleration vector and M, K, and F are the global mass matrix, stiffness matrix, and load matrix, respectively.A direct time-integration scheme is used for solving (4), in which the acceleration and the velocity within each time step are taken as the following equations of the Newmark assumption [18]:(5)D¨n+1=1αΔt2Dn+1-Dn-1αΔtD˙n-12α-1D¨n,D˙n+1=δαΔtDn+1-Dn+1-δαD˙n+1-δ2αΔtD¨n.While denoting V0=D˙(0), A0=D¨(0), and ΔD=Dn+1-Dn, it can be regarded that V0, A0, and ΔD are the initial velocity vector, the initial acceleration vector, and the displacement increment vector of the block system at the beginning of the current time step, respectively. The kinetic energy dissipation of DDA block system depends on the value of gg  (gg=0~1), which is the dynamic controlling parameter. Therefore, gg is introduced to the Newmark assumption; we have(6)D¨n+1=1αΔt2ΔD-gg1αΔtV0-12α-1A0,(7)D˙n+1=δαΔtΔD+gg1-δαV0+1-δ2αΔtA0.Substituting (6) into (4), the motion equation of the block system of the current time step can be derived as(8)K~ΔD=F~,where (9)K~=K+MαΔt2,F~=F+Mgg1αΔtV0+12α-1A0.


## 2.3. Algorithmic Damping of the Newmark Method

For analysis of the stability and accuracy of Newmark method, takingF=0, the following equation can be obtained based on (4)-(5): (10)1+αΩ2Dn+1=2-12-2α+δΩ2Dn-1+12+α-δΩ2Dn-1,where sampling frequency Ω=Δtw and undamped frequency of vibration w=K/M.The solutionsλ of the characteristic equation about (10) are(11)λ1,2=2-g±2-g2-41+h2,where (12)g=1/2+δΩ21+αΩ2,h=1/2-δΩ21+αΩ2.The real solution must have oscillation in the case of small damping, and the stability must be not infinitely growing; thus(13)41+h>2-g2,λ=1+h≤1when Δt is not restricted, (13) should be met, and then α≥(1/4)1/2+δ2, δ≥1/2, 1/2-δ+α≥0.So the conditions of the unconditionally stable Newmark method are(14)δ≥12,α≥1412+δ2.It is thus obvious that [13, 14, 16, 18] both the average acceleration method (α=1/4, δ=1/2) and the constant acceleration method (α=1/2, δ=1) meet the unconditionally stable condition, while the linear acceleration method (α=1/6, δ=1/2) and the central difference method (α=0, δ=1/2) are the conditionally stable method.In general, the accuracy of Newmark method depends on the time interval, the physical parameters of the system, and the loading condition. Here, we only analyze some widely used unconditionally stable integration schemes, to obtain an idea of the accuracy of these methods evaluated based on (11) λ=1+h. Figure 1 shows the algorithmic damping of the Newmark method versus sampling frequency Ω.Figure 1
Accuracy analysis of two Newmark integration schemes.As is shown in Figure1, λ of the average acceleration method is unity, which states that this method has no algorithmic damping. By contrast, λ of the constant acceleration method decreases with the increasing of sampling frequency Ω, which indicates that the constant acceleration method has a rapid attenuation under high frequency interference, while it has very little effect on low frequency. In fact, time interval used in Newmark method is greater than the period of the highest natural frequency of the system; thus, the frequency response is not reliable. So the algorithmic damping is advantageous under certain condition; it may be more suitable for simulating the engineering problems.


## 3. Two-Block Contact Model and the Submatrices

The DDA method is initially developed for discrete block system. For the fundamental study on the damping of DDA, two-block contact model is established, which involves the contact and aims to eliminate the influence of some numeric control parameters (e.g.,g2).


### 3.1. Contact Constraint Conditions of DDA

The DDA method can simulate large displacements of blocks and allows blocks to slide, rotate, and open. A key of the DDA method is contact problem; within each time step, no penetration and no tension between the block interfaces must be satisfied. For 2D DDA method, the contact types include vertex-to-vertex, vertex-to-edge, and edge-to-edge contacts. The edge-to-edge contact can usually be converted to two vertex-to-edge contacts in DDA original program. In this paper, the contact of two blocks is locked, and a normal spring and a shear spring are applied on the contact interface of the two blocks.


### 3.2. Two-Block Contact Model

The model of two contact blocks loaded with the normal force applied along one side of block boundary as shown in Figure2 is established to analyze the damping problem quantitatively. In the model, block B is totally fixed, and block A is subjected to the uniform distributed load on the upper boundary. The length of block A is a and the gravity as well as the contact friction is not considered.Figure 2
Analytical model for DDA damping.Based on the Newmark method and the theory of DDA, the motion equations of two-block system can be established, and then the relationship of some parameters and the influence of damping can be obtained.


### 3.3. Calculations of the Submatrices

As block B is totally fixed, the deformation matrix of block B should be taken as0. The deformation matrix of block A denotes ΔDA=ΔuAΔvAΔrAΔεAxΔεAyΔγAxy. It is clear that (ΔuA,ΔrA,ΔγAxy)=(0,0,0), and the initial deformation matrix of block A is 0. So the global mass matrix, stiffness matrix, and load matrix of two-block system are given as M=MA, K~=KAA, and F~=FA.In the performance of DDA calculation, several submatrices, such as the inertia force submatrix, elastic submatrix, spring stiffness submatrix, load submatrix, and initial stress submatrix, should be deduced and assembled.


#### 3.3.1. Inertia Force Submatrix

The mass matrix of block A is(15)MA=∬ρTATTAdxdy,where ρ is the unit mass.(16)∬ρTATTAdxdy=s000000s000000s1+s2-s3s3s2-s2200-s3s10s3200s30s2s3200s2-s22s32s32s2+s24,where s is the area of block A. s1=sxx-x0sx, s2=syy-y0sy, and s3=sxy-x0sy, where sx is the integral of x over block, sy is the integral of y over block, sxx is the integral of x2 over block, syy is the integral of y2 over block, and sxy is the integral of xy over block.For block A,s=a2, sx=(1/2)a3, sy=(1/2)a3, sxx=(1/3)a4, syy=(1/3)a4, sxy=(1/4)a4, s1=(1/12)a4, s2=(1/12)a4, and s3=0. Therefore, (15) can be rewritten as(17)∬ρTATTAdxdy=ρa2000000a2000000a46000000a412000000a412000000a424.The corresponding inertia force submatrix can be obtained as(18)MαΔt2=1αΔt2∬ρTATTAdxdy⟶KAA,(19)Mgg1αΔtV0+12α-1A0⟶FA.


#### 3.3.2. Elastic Submatrix

The strain energy∏e of block A is(20)∏e=∬12εxεyεxyσxσyτxydxdy=s2∬DATEADAdxdy,where EA is the expanded elastic matrix.The corresponding elastic submatrix can be obtained by minimizing the strain energy∏e(21)sEA=0000000000000000000001μ0000μ10000001-μ2⟶KAA,where E and μ are the elastic modulus and the poisson ratio, respectively.


#### 3.3.3. Spring Stiffness Submatrix

The stiffness submatrix of normal spring [2] is(22)per,Mer,MT⟶KAA,-ps0er,Ml⟶FA,where l=a, er,M=[(y2-y3)t1r(x1,y1)+(x3-x2)t2r(x1,y1)]/l, and (23)s0=1x1y11x2y21x3y3.As shown in Figure 2, (x1,y1), (x2,y2), and (x3,y3) are (0,0), (0,0), and (a,0), respectively. Consequently er,M is simplified as(24)er,M=01-a20-a2-a4T.Similarly, for the contact point N, (x1,y1), (x2,y2), and (x3,y3) are (a,0), (0,0), and (a,0), respectively. Thus er,N of block A is er,N=01a/20-a/2a/4T.The stiffness submatrix [2] of shear spring is(25)per,M′er,M′T⟶KAA,-ps0er,M′l⟶FA,where(26)er,M′=x3-x2t1rx1,y1+y3-y2t2rx1,y1l.Similarly, er,M′=10a/2-a/20-a/4T, er,N′=10a/2a/20-a/4T.As a result, the spring stiffness submatrix of blockA can be acquired as(27)per,Mer,MT+per,Ner,NT+per,M′er,M′T+per,N′er,N′T=p20a00-a20200-a0a0a2000000a22000-a00a220-a20000a24⟶KAA,(28)-ps0er,Ml-ps0er,Nl-ps0er,M′l-ps0er,N′l⟶FA.


#### 3.3.4. Load Submatrix

As is shown in Figure2, the uniform load Fx,Fy is equal to 0,q. The corresponding load submatrix [2] of block A can be presented as(29)∫01TiTFxFyldt=lFxFy-12y2+y1-2y0Fx+12x2+x1-2x0Fy12x2+x1-2x0Fx12y2+y1-2y0Fy14y2+y1-2y0Fx+14x2+x1-2x0Fy=0-aq00-a2q20⟶FA.


#### 3.3.5. Initial Stress Submatrix

The corresponding initial stress submatrix is(30)-s000σAx,nσAy,nτAxy,n⟶FA,where σAx,n, σAy,n, and τAxy,n are the stress components of nth time step of block A.


## 3.1. Contact Constraint Conditions of DDA

The DDA method can simulate large displacements of blocks and allows blocks to slide, rotate, and open. A key of the DDA method is contact problem; within each time step, no penetration and no tension between the block interfaces must be satisfied. For 2D DDA method, the contact types include vertex-to-vertex, vertex-to-edge, and edge-to-edge contacts. The edge-to-edge contact can usually be converted to two vertex-to-edge contacts in DDA original program. In this paper, the contact of two blocks is locked, and a normal spring and a shear spring are applied on the contact interface of the two blocks.


## 3.2. Two-Block Contact Model

The model of two contact blocks loaded with the normal force applied along one side of block boundary as shown in Figure2 is established to analyze the damping problem quantitatively. In the model, block B is totally fixed, and block A is subjected to the uniform distributed load on the upper boundary. The length of block A is a and the gravity as well as the contact friction is not considered.Figure 2
Analytical model for DDA damping.Based on the Newmark method and the theory of DDA, the motion equations of two-block system can be established, and then the relationship of some parameters and the influence of damping can be obtained.


## 3.3. Calculations of the Submatrices

As block B is totally fixed, the deformation matrix of block B should be taken as0. The deformation matrix of block A denotes ΔDA=ΔuAΔvAΔrAΔεAxΔεAyΔγAxy. It is clear that (ΔuA,ΔrA,ΔγAxy)=(0,0,0), and the initial deformation matrix of block A is 0. So the global mass matrix, stiffness matrix, and load matrix of two-block system are given as M=MA, K~=KAA, and F~=FA.In the performance of DDA calculation, several submatrices, such as the inertia force submatrix, elastic submatrix, spring stiffness submatrix, load submatrix, and initial stress submatrix, should be deduced and assembled.


### 3.3.1. Inertia Force Submatrix

The mass matrix of block A is(15)MA=∬ρTATTAdxdy,where ρ is the unit mass.(16)∬ρTATTAdxdy=s000000s000000s1+s2-s3s3s2-s2200-s3s10s3200s30s2s3200s2-s22s32s32s2+s24,where s is the area of block A. s1=sxx-x0sx, s2=syy-y0sy, and s3=sxy-x0sy, where sx is the integral of x over block, sy is the integral of y over block, sxx is the integral of x2 over block, syy is the integral of y2 over block, and sxy is the integral of xy over block.For block A,s=a2, sx=(1/2)a3, sy=(1/2)a3, sxx=(1/3)a4, syy=(1/3)a4, sxy=(1/4)a4, s1=(1/12)a4, s2=(1/12)a4, and s3=0. Therefore, (15) can be rewritten as(17)∬ρTATTAdxdy=ρa2000000a2000000a46000000a412000000a412000000a424.The corresponding inertia force submatrix can be obtained as(18)MαΔt2=1αΔt2∬ρTATTAdxdy⟶KAA,(19)Mgg1αΔtV0+12α-1A0⟶FA.


### 3.3.2. Elastic Submatrix

The strain energy∏e of block A is(20)∏e=∬12εxεyεxyσxσyτxydxdy=s2∬DATEADAdxdy,where EA is the expanded elastic matrix.The corresponding elastic submatrix can be obtained by minimizing the strain energy∏e(21)sEA=0000000000000000000001μ0000μ10000001-μ2⟶KAA,where E and μ are the elastic modulus and the poisson ratio, respectively.


### 3.3.3. Spring Stiffness Submatrix

The stiffness submatrix of normal spring [2] is(22)per,Mer,MT⟶KAA,-ps0er,Ml⟶FA,where l=a, er,M=[(y2-y3)t1r(x1,y1)+(x3-x2)t2r(x1,y1)]/l, and (23)s0=1x1y11x2y21x3y3.As shown in Figure 2, (x1,y1), (x2,y2), and (x3,y3) are (0,0), (0,0), and (a,0), respectively. Consequently er,M is simplified as(24)er,M=01-a20-a2-a4T.Similarly, for the contact point N, (x1,y1), (x2,y2), and (x3,y3) are (a,0), (0,0), and (a,0), respectively. Thus er,N of block A is er,N=01a/20-a/2a/4T.The stiffness submatrix [2] of shear spring is(25)per,M′er,M′T⟶KAA,-ps0er,M′l⟶FA,where(26)er,M′=x3-x2t1rx1,y1+y3-y2t2rx1,y1l.Similarly, er,M′=10a/2-a/20-a/4T, er,N′=10a/2a/20-a/4T.As a result, the spring stiffness submatrix of blockA can be acquired as(27)per,Mer,MT+per,Ner,NT+per,M′er,M′T+per,N′er,N′T=p20a00-a20200-a0a0a2000000a22000-a00a220-a20000a24⟶KAA,(28)-ps0er,Ml-ps0er,Nl-ps0er,M′l-ps0er,N′l⟶FA.


### 3.3.4. Load Submatrix

As is shown in Figure2, the uniform load Fx,Fy is equal to 0,q. The corresponding load submatrix [2] of block A can be presented as(29)∫01TiTFxFyldt=lFxFy-12y2+y1-2y0Fx+12x2+x1-2x0Fy12x2+x1-2x0Fx12y2+y1-2y0Fy14y2+y1-2y0Fx+14x2+x1-2x0Fy=0-aq00-a2q20⟶FA.


### 3.3.5. Initial Stress Submatrix

The corresponding initial stress submatrix is(30)-s000σAx,nσAy,nτAxy,n⟶FA,where σAx,n, σAy,n, and τAxy,n are the stress components of nth time step of block A.


## 3.3.1. Inertia Force Submatrix

The mass matrix of block A is(15)MA=∬ρTATTAdxdy,where ρ is the unit mass.(16)∬ρTATTAdxdy=s000000s000000s1+s2-s3s3s2-s2200-s3s10s3200s30s2s3200s2-s22s32s32s2+s24,where s is the area of block A. s1=sxx-x0sx, s2=syy-y0sy, and s3=sxy-x0sy, where sx is the integral of x over block, sy is the integral of y over block, sxx is the integral of x2 over block, syy is the integral of y2 over block, and sxy is the integral of xy over block.For block A,s=a2, sx=(1/2)a3, sy=(1/2)a3, sxx=(1/3)a4, syy=(1/3)a4, sxy=(1/4)a4, s1=(1/12)a4, s2=(1/12)a4, and s3=0. Therefore, (15) can be rewritten as(17)∬ρTATTAdxdy=ρa2000000a2000000a46000000a412000000a412000000a424.The corresponding inertia force submatrix can be obtained as(18)MαΔt2=1αΔt2∬ρTATTAdxdy⟶KAA,(19)Mgg1αΔtV0+12α-1A0⟶FA.


## 3.3.2. Elastic Submatrix

The strain energy∏e of block A is(20)∏e=∬12εxεyεxyσxσyτxydxdy=s2∬DATEADAdxdy,where EA is the expanded elastic matrix.The corresponding elastic submatrix can be obtained by minimizing the strain energy∏e(21)sEA=0000000000000000000001μ0000μ10000001-μ2⟶KAA,where E and μ are the elastic modulus and the poisson ratio, respectively.


## 3.3.3. Spring Stiffness Submatrix

The stiffness submatrix of normal spring [2] is(22)per,Mer,MT⟶KAA,-ps0er,Ml⟶FA,where l=a, er,M=[(y2-y3)t1r(x1,y1)+(x3-x2)t2r(x1,y1)]/l, and (23)s0=1x1y11x2y21x3y3.As shown in Figure 2, (x1,y1), (x2,y2), and (x3,y3) are (0,0), (0,0), and (a,0), respectively. Consequently er,M is simplified as(24)er,M=01-a20-a2-a4T.Similarly, for the contact point N, (x1,y1), (x2,y2), and (x3,y3) are (a,0), (0,0), and (a,0), respectively. Thus er,N of block A is er,N=01a/20-a/2a/4T.The stiffness submatrix [2] of shear spring is(25)per,M′er,M′T⟶KAA,-ps0er,M′l⟶FA,where(26)er,M′=x3-x2t1rx1,y1+y3-y2t2rx1,y1l.Similarly, er,M′=10a/2-a/20-a/4T, er,N′=10a/2a/20-a/4T.As a result, the spring stiffness submatrix of blockA can be acquired as(27)per,Mer,MT+per,Ner,NT+per,M′er,M′T+per,N′er,N′T=p20a00-a20200-a0a0a2000000a22000-a00a220-a20000a24⟶KAA,(28)-ps0er,Ml-ps0er,Nl-ps0er,M′l-ps0er,N′l⟶FA.


## 3.3.4. Load Submatrix

As is shown in Figure2, the uniform load Fx,Fy is equal to 0,q. The corresponding load submatrix [2] of block A can be presented as(29)∫01TiTFxFyldt=lFxFy-12y2+y1-2y0Fx+12x2+x1-2x0Fy12x2+x1-2x0Fx12y2+y1-2y0Fy14y2+y1-2y0Fx+14x2+x1-2x0Fy=0-aq00-a2q20⟶FA.


## 3.3.5. Initial Stress Submatrix

The corresponding initial stress submatrix is(30)-s000σAx,nσAy,nτAxy,n⟶FA,where σAx,n, σAy,n, and τAxy,n are the stress components of nth time step of block A.


## 4. Damping Analysis of DDA

Based on the two-block contact model, the influencing factors of damping can be analyzed quantitatively. By assembling the global stiffness matrix and the load matrix, the total equation motion can be obtained.


### 4.1. Formula of the Total Equation Motion

Solving (8) based on two-block model, K~ is formed by (18), (21), and (27), F~ is formed by (19) and (28)–(30). As (ΔuA,ΔrA,ΔγAxy)=(0,0,0), the following three iteration equations can be obtained:(31)ρa2αΔt2+2pΔvA,n+1+-paΔεAy,n+1=f1,n+1,ρa412αΔt2+Ea21-μ2+a2p2ΔεAx,n+1+Eμa21-μ2ΔεAy,n+1=f2,n+1,-paΔvA,n+1+Eμa21-μ2ΔεAx,n+1+ρa412αΔt2+Ea21-μ2+a2p2ΔεAy,n+1=f3,n+1,where(32)f1,n+1=-pvM,n+vN,n-aq+ggρa2αΔtV1,n+ρa212α-1A1,n,f2,n+1=-a2σAx,n+ggρa412αΔtV2,n+ρa41212α-1A2,n,f3,n+1=pa2vM,n+vN,n-a2σAy,n-a2q2+ggρa412αΔtV3,n+ρa41212α-1A3,n,ΔvA,n+1, ΔεAx,n+1, and ΔεAy,n+1 are the deformation increments at the n+1th time step of block A, V0=0V1,n0V2,nV3,n0 and A0=0A1,n0A2,nA3,n0 are the initial velocity and initial acceleration at the n+1th time step of block A, vM,n and vN,n are the accumulated displacement at points M and N, and σAx,n and σAy,n are the stress components at the nth time step.After the valuesα, δ, and gg are determined, using the initial conditions, (V1,1,V2,1,V3,1)=(0,0,0), (A1,1,A2,1,A3,1)=(0,0,0), (σAx,1,σAy,1)=(0,0), and (vM,1,vN,1)=(0,0), ΔvA,n+1, ΔεAx,n+1, and ΔεAy,n+1 can be obtained by solving (31). Substituting (ΔvA,n+1,ΔεAx,n+1,ΔεAy,n+1) into (2), the physical equations of elastic dynamics, (6) and (7), respectively, (vM,n+1,vN,n+1), (σAx,n+1,σAy,n+1), (A1,n+1,A2,n+1,A3,n+1), and (A1,n+1,A2,n+1,A3,n+1) can be obtained.


### 4.2. Algorithmic Damping of DDA

#### 4.2.1. Influencing Factors of Algorithmic Damping

For studying the algorithmic damping and eliminating the influence of kinetic damping, the kinetic damping is not considered by denotinggg=1. The acceleration within each time step of DDA is assumed to be constant (α=1/2, δ=1).Considering thatρ, E, μ, and q are the physical and mechanical parameters, we only study the influences of time interval Δt and the spring stiffness p. In order to perform case studies, assumptions are made that a=2.0 m, ρ=2700 kg/m3, E=1 GPa, μ=0, and q=5 MPa. The curve of displacement of the middle point of block A versus time step under the different time interval Δt (0.01 sec, 0.002 sec, 0.001 sec) is shown in Figure 3. The corresponding contact force curve is presented in Figure 4. The curve of displacement of the middle point of block A versus time step under the different spring stiffness p (1E, 10E, 100E) is shown in Figure 5. The corresponding contact force curve is shown in Figure 6. The results show that contact force approximates the real solution 5 MPa successively and keeps stable at last, which satisfies the equilibrium of force system.Figure 3
Displacement of the middle point versus time steps(p=1E).Figure 4
Contact force versus time steps(p=1E).Figure 5
Displacement of the middle point versus time steps(Δt=0.002 sec).Figure 6
Contact force versus time steps(Δt=0.002 sec).As shown in Figures3 and 4, the smaller the time interval is, the more powerful the displacement and contact force oscillate. With the increasing of time interval, the displacement and contact force can converge to the equilibrium of position rapidly, which is actually caused by the algorithmic damping. For example, if Δt=0.01 sec, the displacement and contact force reach their convergence state after five and four time steps, respectively. But if Δt=0.002 sec, the time consumption for convergence of displacement and contact force is thirty-four time steps and thirty-six time steps, respectively, much more than those of Δt=0.01 sec.Clearly, from Figures5 and 6, as the value of spring stiffness decreases, the oscillation becomes more obvious. However, when spring stiffness increases to a certain degree, it causes no difference. For example, if p=1E, the convergence of both displacement and contact force costs forty-four time steps. By contrast, if p=10E, the convergence of displacement and contact force costs twenty and twenty-one time steps, respectively. In addition, the difference in convergence of p=100E and p=10E is not obvious.Furthermore, when time interval is 0.00001 sec, the relationship between displacement of the middle point and time steps tends to be vibration curves with no attenuation, and spring stiffness only influences the amplitude and frequency of vibration curves, which can be seen in Figure7.Figure 7
Displacement of the middle point versus time steps(Δt=0.00001 sec).


#### 4.2.2. Analysis of Algorithmic Damping

The deformation and failure of rocks is a process of energy dissipation; the damping research of DDA is a fundamental problem. Algorithm damping can be caused by DDA using Newmark constant acceleration integral scheme. If average acceleration method is used in DDA simulation, there will be no algorithmic damping, which can be seen in Figure8. Significantly, the algorithm damping is useful under certain conditions, it can guarantee the open-close convergence, and for the quasi-static problem, it can quickly converge to the correct solution.Figure 8
The relative error versus time steps(Δt=0.001 sec, p=1E).From Section4.2.1, it can be seen that time interval is a critical factor of the algorithmic damping. It has been shown that, with increasing the time interval, the displacement and contact force can converge to the equilibrium of position rapidly. However, time interval is an artificial parameter because the values have significant effect on the accuracy of the results. For instance, if the time interval is too large, it may cause the loss of diagonal dominance, interfering with the accuracy. If the time interval is too small, the algorithmic damping may be neglected; thus the contact displacement may vibrate indefinitely, which seriously affects the convergence of open-close iteration.Moreover, the smaller the value of spring stiffness is, the more obvious the oscillation is, and when spring stiffness is large enough, it causes no difference. These phenomena illustrate that the algorithmic damping decreases with the decreasing of spring stiffness under certain conditions. Additionally, when the time interval is too small, the spring stiffness may have no obvious effect on the algorithmic damping and only influence the amplitude and frequency of vibration curve, which explains that time interval is more sensitive to the algorithmic damping.


### 4.3. Kinetic Damping of DDA

For a better understanding of the kinetic damping, the acceleration within each time step of DDA is assumed to be average acceleration whereα=1/4, δ=1/2; in that case the algorithmic damping of DDA can be neglected. The physical and mechanical parameters are the same as Section 4.2.1. The curve of displacement of the middle point of block A versus time steps under the different gg (1.00, 0.98, 0.95, 0.90) is shown in Figure 9.Figure 9
Displacement amplitude attenuation of different values ofgg(Δt=0.001 sec, p=1E).Some scholars have researched the kinetic damping, and most of them made their efforts to compare the simulation results with physical-experimental results, which aimed to obtain the discounted values of the initial velocity of each time step. However, because of the complexity of the energy dissipation, the accurate discounted value is still difficult to determine. Here, we investigate the relative trend of kinetic damping under the different values of dynamic control parameter, where the corresponding discounted values are 0%, 2%, 5%, and 1%. TakingΔt=0.001 sec and p=1E, curves of the displacement against time steps can be obtained as given in Figure 9, from which the amplitude of vibration curve is found to attenuate fast with the increase of the discounted value.


## 4.1. Formula of the Total Equation Motion

Solving (8) based on two-block model, K~ is formed by (18), (21), and (27), F~ is formed by (19) and (28)–(30). As (ΔuA,ΔrA,ΔγAxy)=(0,0,0), the following three iteration equations can be obtained:(31)ρa2αΔt2+2pΔvA,n+1+-paΔεAy,n+1=f1,n+1,ρa412αΔt2+Ea21-μ2+a2p2ΔεAx,n+1+Eμa21-μ2ΔεAy,n+1=f2,n+1,-paΔvA,n+1+Eμa21-μ2ΔεAx,n+1+ρa412αΔt2+Ea21-μ2+a2p2ΔεAy,n+1=f3,n+1,where(32)f1,n+1=-pvM,n+vN,n-aq+ggρa2αΔtV1,n+ρa212α-1A1,n,f2,n+1=-a2σAx,n+ggρa412αΔtV2,n+ρa41212α-1A2,n,f3,n+1=pa2vM,n+vN,n-a2σAy,n-a2q2+ggρa412αΔtV3,n+ρa41212α-1A3,n,ΔvA,n+1, ΔεAx,n+1, and ΔεAy,n+1 are the deformation increments at the n+1th time step of block A, V0=0V1,n0V2,nV3,n0 and A0=0A1,n0A2,nA3,n0 are the initial velocity and initial acceleration at the n+1th time step of block A, vM,n and vN,n are the accumulated displacement at points M and N, and σAx,n and σAy,n are the stress components at the nth time step.After the valuesα, δ, and gg are determined, using the initial conditions, (V1,1,V2,1,V3,1)=(0,0,0), (A1,1,A2,1,A3,1)=(0,0,0), (σAx,1,σAy,1)=(0,0), and (vM,1,vN,1)=(0,0), ΔvA,n+1, ΔεAx,n+1, and ΔεAy,n+1 can be obtained by solving (31). Substituting (ΔvA,n+1,ΔεAx,n+1,ΔεAy,n+1) into (2), the physical equations of elastic dynamics, (6) and (7), respectively, (vM,n+1,vN,n+1), (σAx,n+1,σAy,n+1), (A1,n+1,A2,n+1,A3,n+1), and (A1,n+1,A2,n+1,A3,n+1) can be obtained.


## 4.2. Algorithmic Damping of DDA

### 4.2.1. Influencing Factors of Algorithmic Damping

For studying the algorithmic damping and eliminating the influence of kinetic damping, the kinetic damping is not considered by denotinggg=1. The acceleration within each time step of DDA is assumed to be constant (α=1/2, δ=1).Considering thatρ, E, μ, and q are the physical and mechanical parameters, we only study the influences of time interval Δt and the spring stiffness p. In order to perform case studies, assumptions are made that a=2.0 m, ρ=2700 kg/m3, E=1 GPa, μ=0, and q=5 MPa. The curve of displacement of the middle point of block A versus time step under the different time interval Δt (0.01 sec, 0.002 sec, 0.001 sec) is shown in Figure 3. The corresponding contact force curve is presented in Figure 4. The curve of displacement of the middle point of block A versus time step under the different spring stiffness p (1E, 10E, 100E) is shown in Figure 5. The corresponding contact force curve is shown in Figure 6. The results show that contact force approximates the real solution 5 MPa successively and keeps stable at last, which satisfies the equilibrium of force system.Figure 3
Displacement of the middle point versus time steps(p=1E).Figure 4
Contact force versus time steps(p=1E).Figure 5
Displacement of the middle point versus time steps(Δt=0.002 sec).Figure 6
Contact force versus time steps(Δt=0.002 sec).As shown in Figures3 and 4, the smaller the time interval is, the more powerful the displacement and contact force oscillate. With the increasing of time interval, the displacement and contact force can converge to the equilibrium of position rapidly, which is actually caused by the algorithmic damping. For example, if Δt=0.01 sec, the displacement and contact force reach their convergence state after five and four time steps, respectively. But if Δt=0.002 sec, the time consumption for convergence of displacement and contact force is thirty-four time steps and thirty-six time steps, respectively, much more than those of Δt=0.01 sec.Clearly, from Figures5 and 6, as the value of spring stiffness decreases, the oscillation becomes more obvious. However, when spring stiffness increases to a certain degree, it causes no difference. For example, if p=1E, the convergence of both displacement and contact force costs forty-four time steps. By contrast, if p=10E, the convergence of displacement and contact force costs twenty and twenty-one time steps, respectively. In addition, the difference in convergence of p=100E and p=10E is not obvious.Furthermore, when time interval is 0.00001 sec, the relationship between displacement of the middle point and time steps tends to be vibration curves with no attenuation, and spring stiffness only influences the amplitude and frequency of vibration curves, which can be seen in Figure7.Figure 7
Displacement of the middle point versus time steps(Δt=0.00001 sec).


### 4.2.2. Analysis of Algorithmic Damping

The deformation and failure of rocks is a process of energy dissipation; the damping research of DDA is a fundamental problem. Algorithm damping can be caused by DDA using Newmark constant acceleration integral scheme. If average acceleration method is used in DDA simulation, there will be no algorithmic damping, which can be seen in Figure8. Significantly, the algorithm damping is useful under certain conditions, it can guarantee the open-close convergence, and for the quasi-static problem, it can quickly converge to the correct solution.Figure 8
The relative error versus time steps(Δt=0.001 sec, p=1E).From Section4.2.1, it can be seen that time interval is a critical factor of the algorithmic damping. It has been shown that, with increasing the time interval, the displacement and contact force can converge to the equilibrium of position rapidly. However, time interval is an artificial parameter because the values have significant effect on the accuracy of the results. For instance, if the time interval is too large, it may cause the loss of diagonal dominance, interfering with the accuracy. If the time interval is too small, the algorithmic damping may be neglected; thus the contact displacement may vibrate indefinitely, which seriously affects the convergence of open-close iteration.Moreover, the smaller the value of spring stiffness is, the more obvious the oscillation is, and when spring stiffness is large enough, it causes no difference. These phenomena illustrate that the algorithmic damping decreases with the decreasing of spring stiffness under certain conditions. Additionally, when the time interval is too small, the spring stiffness may have no obvious effect on the algorithmic damping and only influence the amplitude and frequency of vibration curve, which explains that time interval is more sensitive to the algorithmic damping.


## 4.2.1. Influencing Factors of Algorithmic Damping

For studying the algorithmic damping and eliminating the influence of kinetic damping, the kinetic damping is not considered by denotinggg=1. The acceleration within each time step of DDA is assumed to be constant (α=1/2, δ=1).Considering thatρ, E, μ, and q are the physical and mechanical parameters, we only study the influences of time interval Δt and the spring stiffness p. In order to perform case studies, assumptions are made that a=2.0 m, ρ=2700 kg/m3, E=1 GPa, μ=0, and q=5 MPa. The curve of displacement of the middle point of block A versus time step under the different time interval Δt (0.01 sec, 0.002 sec, 0.001 sec) is shown in Figure 3. The corresponding contact force curve is presented in Figure 4. The curve of displacement of the middle point of block A versus time step under the different spring stiffness p (1E, 10E, 100E) is shown in Figure 5. The corresponding contact force curve is shown in Figure 6. The results show that contact force approximates the real solution 5 MPa successively and keeps stable at last, which satisfies the equilibrium of force system.Figure 3
Displacement of the middle point versus time steps(p=1E).Figure 4
Contact force versus time steps(p=1E).Figure 5
Displacement of the middle point versus time steps(Δt=0.002 sec).Figure 6
Contact force versus time steps(Δt=0.002 sec).As shown in Figures3 and 4, the smaller the time interval is, the more powerful the displacement and contact force oscillate. With the increasing of time interval, the displacement and contact force can converge to the equilibrium of position rapidly, which is actually caused by the algorithmic damping. For example, if Δt=0.01 sec, the displacement and contact force reach their convergence state after five and four time steps, respectively. But if Δt=0.002 sec, the time consumption for convergence of displacement and contact force is thirty-four time steps and thirty-six time steps, respectively, much more than those of Δt=0.01 sec.Clearly, from Figures5 and 6, as the value of spring stiffness decreases, the oscillation becomes more obvious. However, when spring stiffness increases to a certain degree, it causes no difference. For example, if p=1E, the convergence of both displacement and contact force costs forty-four time steps. By contrast, if p=10E, the convergence of displacement and contact force costs twenty and twenty-one time steps, respectively. In addition, the difference in convergence of p=100E and p=10E is not obvious.Furthermore, when time interval is 0.00001 sec, the relationship between displacement of the middle point and time steps tends to be vibration curves with no attenuation, and spring stiffness only influences the amplitude and frequency of vibration curves, which can be seen in Figure7.Figure 7
Displacement of the middle point versus time steps(Δt=0.00001 sec).


## 4.2.2. Analysis of Algorithmic Damping

The deformation and failure of rocks is a process of energy dissipation; the damping research of DDA is a fundamental problem. Algorithm damping can be caused by DDA using Newmark constant acceleration integral scheme. If average acceleration method is used in DDA simulation, there will be no algorithmic damping, which can be seen in Figure8. Significantly, the algorithm damping is useful under certain conditions, it can guarantee the open-close convergence, and for the quasi-static problem, it can quickly converge to the correct solution.Figure 8
The relative error versus time steps(Δt=0.001 sec, p=1E).From Section4.2.1, it can be seen that time interval is a critical factor of the algorithmic damping. It has been shown that, with increasing the time interval, the displacement and contact force can converge to the equilibrium of position rapidly. However, time interval is an artificial parameter because the values have significant effect on the accuracy of the results. For instance, if the time interval is too large, it may cause the loss of diagonal dominance, interfering with the accuracy. If the time interval is too small, the algorithmic damping may be neglected; thus the contact displacement may vibrate indefinitely, which seriously affects the convergence of open-close iteration.Moreover, the smaller the value of spring stiffness is, the more obvious the oscillation is, and when spring stiffness is large enough, it causes no difference. These phenomena illustrate that the algorithmic damping decreases with the decreasing of spring stiffness under certain conditions. Additionally, when the time interval is too small, the spring stiffness may have no obvious effect on the algorithmic damping and only influence the amplitude and frequency of vibration curve, which explains that time interval is more sensitive to the algorithmic damping.


## 4.3. Kinetic Damping of DDA

For a better understanding of the kinetic damping, the acceleration within each time step of DDA is assumed to be average acceleration whereα=1/4, δ=1/2; in that case the algorithmic damping of DDA can be neglected. The physical and mechanical parameters are the same as Section 4.2.1. The curve of displacement of the middle point of block A versus time steps under the different gg (1.00, 0.98, 0.95, 0.90) is shown in Figure 9.Figure 9
Displacement amplitude attenuation of different values ofgg(Δt=0.001 sec, p=1E).Some scholars have researched the kinetic damping, and most of them made their efforts to compare the simulation results with physical-experimental results, which aimed to obtain the discounted values of the initial velocity of each time step. However, because of the complexity of the energy dissipation, the accurate discounted value is still difficult to determine. Here, we investigate the relative trend of kinetic damping under the different values of dynamic control parameter, where the corresponding discounted values are 0%, 2%, 5%, and 1%. TakingΔt=0.001 sec and p=1E, curves of the displacement against time steps can be obtained as given in Figure 9, from which the amplitude of vibration curve is found to attenuate fast with the increase of the discounted value.


## 5. Conclusions

The analysis of the seismic response is of great significance in rock engineering. The DDA method is very suitable for the problem, which has been confirmed by many scholars. However, the deformation and failure of rocks is a process of energy dissipation; the damping research of DDA is a fundamental problem. In this paper, we perform quantitative analysis of the influencing factors of damping based on the Newmark method and the mechanical model of two contact blocks. Through the analysis of damping, the following understandings are obtained.(1) The time interval is a controlling factor of the algorithmic damping. The algorithmic damping increases with the increasing of time interval; the displacement and contact force can converge to the equilibrium of position rapidly with larger time interval.(2) The spring stiffness has a certain effect on the algorithmic damping. The algorithmic damping increases as spring stiffness increases, but if spring stiffness is large enough, it will cause no difference.(3) Given a small enough time interval, the algorithmic damping may be neglected, and the spring stiffness may have no obvious effect on the algorithmic damping; it can only influence the amplitude and frequency of vibration curve. It is clear that the time interval is more sensitive to algorithmic damping. The essence of time interval influencing the open-close iteration is the algorithmic damping mainly controlled by the time interval.(4) The kinetic damping of DDA is only dependent on the discounted values of the initial velocity of each time step. The amplitude of vibration curve attenuates fast with the increase of the discounted value, which illustrates that the appropriate discounted value is critical to the energy dissipation.



---

*Source: 7480937-2017-07-18.xml*