# Improved DDA Method Based on Explicitly Solving Contact Constraints

Authors: Jian Zhao; Ming Xiao; Juntao Chen; Dongdong Li
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This paper proposes an improved DDA method based on explicitly solving contact constraints. The potential energy function generated by contacting, which contains only displacement variables as an unknown, is deduced based on the approximated step function and Lagrange interpolation, and the displacement variables and contact constraints are obtained via the variable metric method by analyzing the potential energy extremum. There is no need to conduct the open-close iteration during the process of calculation. The improved DDA method based on explicitly solving contact constraints has high precision and a more stable and more robust computational convergence. The accuracy and iterative stability of the improved DDA method are verified using two numerical examples.


---

## Body

## 1. Introduction

The discontinuous deformation analysis method proposed by Shi [1] is a new, efficient numerical method for discontinuous rock mass. The DDA method can efficiently simulate mutual contact between blocks, account for the impact of geological joints on block movement, and calculate the displacement and deformation of rock mass.In the past 20 years, numerous studies have been performed to overcome the limitations of the DDA method. Koo and Max [2, 3] applied high-order polynomials to the displacement functions of blocks to improve the accuracy of the calculation results. Shyu and Yue [4, 5] coupled the finite element method and DDA method and used the expression of the finite element method to describe the fields of displacement and stress in blocks. Amadei et al. [6] divided the original large block into several small blocks and proposed the secondary block theory, which can describe the rupture behavior and its propagation inside blocks. Beyabanaki and Lin [7, 8] used an augmented Lagrange multiplier method instead of the penalty function method to address the contact issue in the DDA method, avoiding the difficulty of selecting the contact stiffness. Koo and Chern [9] proposed the rigid DDA method to reduce the volume expansion error due to rotation. MacLaughlin and Doolin [10] reviewed the validation of the DDA method in the past years. The development of the DDA method has also played a role in other subjects [11].The DDA method also encounters a contact issue. Addressing this issue is important but difficult; whether to consider the contact problem is one of the main features that distinguishes the DDA method from the finite element method [12]. Because the objects analyzed in the DDA method are discontinuous blocks and the number of contact pairs increases significantly with an increasing number of blocks, the problem of mutual contact between blocks is one of the main tasks analyzed during DDA analysis. The penalty function method [13] was introduced in the conventional DDA method; a spring is used to simulate the contact force, which allows for a certain degree of embedding distance between blocks and controls the calculating accuracy through the penalty factor. Amadei et al. [14] introduced the augmented Lagrange multiplier method to realize the contact constraint between adjacent blocks; the algorithm can accurately meet the contact conditions set. The problem that exists in both common methods is that they cannot predict the types of contact constraints between blocks and cannot directly calculate certain contact potential energies. The correct results can be obtained only by modifying the types of contact constraints and the corresponding contact springs repeatedly, which is called open-close iteration. Open-close iteration can be viewed as an implicit algorithm and cannot be used to obtain the types of contact constraints between blocks directly. It is difficult to meet the convergence condition because of the poor convergence of implicit open-close iteration and the stringent time step demand [15]; moreover, implicit open-close iteration is extremely time consuming. When the block model is complex, the types of contact constraints are more difficult to obtain and the time consumed in conducting open-close iteration increases dramatically with an increasing number of contact pairs, which leads to low computational efficiency.Few studies have focused on improving the open-close iteration; instead, the majority of research work is still focused on the study of contact judgment [16–18]. Zheng [19, 20] proposed the discontinuous deformation analysis method based on linear complementary theory, which originates from the variational equations of momentum conservation. The linear complementary relationship is used to describe the normal and tangential contact constraint and avoids the implicit open-close iteration; however, iteration is still needed when determining the valid access line of edge-to-edge contact. Feng and Wu [21, 22] deduced the theoretical formula of the block displacement under the static condition and analyzed the complete implementation process of the open-close iteration in the DDA method, which proved the rationality of the open-close iteration but did not improve the poor convergence of the implicit open-close iteration. Zhang and Cheng [23] proposed the variable stiffness spring method, in which the number of springs is modified to achieve a reasonable stiffness. The separation of shear and normal stiffness for the DDA method was introduced, but it still requires the contact constraint to be modified.This paper proposes an improved DDA method based on explicitly solving contact constraints to improve the convergence in the implicit open-close iteration of discontinuous deformation analysis. The process of conducting the conventional DDA method can be viewed as a process of obtaining the extremum of the potential energy function of the block system with the contact constraints. The contact constraints and displacement variables of blocks are regarded as two groups of unknowns, but they are not completely separated and unrelated. The contact constraints can be represented by displacement variables, and they are interrelated. The relationship between the contact constraints and displacement variables of blocks is studied, and the overall potential energy function of a block system containing only displacement variables as the unknown is deduced, which does not require implicit open-close iteration. The process can be viewed as an explicit algorithm for determining the type of contact constraints and does not require iteration. Then, the minimum extremum of the potential energy function of the block system is obtained using the variable metric method, and the final steady state of the block system is achieved.


## 2. Contact Detection

Each time step has corresponding contact constraints between blocks, and all types of contacts can be grouped as either vertex-edge contacts or vertex-vertex contacts. Cundall [24] proposed the common-plane method in 1988, in which the type of contact can be determined by the number of the points contacting the common plane. All types of contacts can be converted into one or more vertex-edge contacts, and a vertex-vertex contact can be viewed as a vertex-edge contact with an uncertain valid access line. From this perspective, all types of contacts can be seen as vertex-edge contacts. Next, the relationship between embedding and displacement is discussed in detail.Suppose that pointP1 moves to the position of point P1′ after deforming, P2P3 is the contact line, and xi,yi and ui,vi are the coordinate values and displacement of point Pii=1,2,3, respectively. Embedding occurs if points P1, P2, and P3 rotate from 0x to 0y in the same direction and point P1′ crosses over the line P2P3, which is illustrated in Figure 1.Figure 1
Vertex-edge contact in the block.
(a)
Before embedding (b)
After embeddingEmbedding can be expressed as(1)Δ=1x1+u1y1+v11x2+u2y2+v21x3+u3y3+v3<0.Equation (1) illustrates that Δ is a vector-valued function containing only the displacement variables D as an unknown, which can be expressed as (2)Δ=ΔD.IfΔ is smaller than 0, embedding occurs and normal and tangential embedding distances exist. The normal embedding distance dn and the tangential embedding distance ds can be expressed as [1](3)dn=H1TDi+G1TDj+S1l,ds=H2TDi+G2TDj+S2l,where Di and Dj are the displacement variables of the two blocks in contact with each other and H1, H2, G1, G2, S1, S2, and l are the intermediate variables.Equation (3) illustrates that the normal embedding distance dn and the tangential embedding distance ds are vector-valued function containing only the displacement D as an unknown, which can be expressed as(4)dn=dnD,ds=dsD.VariableΔ, normal embedding distance dn, and tangential embedding distance ds are deduced for the following deduction of the approximated step function.


## 3. Contact Constraints Simulated by the Approximated Step Function

All types of contact can be grouped into the following three types: open, sliding, and fixed. The rules of mode changes are shown in Table1.Table 1
Discriminant table of contact constraints.
Time step →i-1 Time step →i Open Sliding Fixed Open d n > 0 d n < 0   a n d   d s > d n tan ⁡ φ d n < 0   a n d   d s < d n tan ⁡ φ Sliding d n > 0 d n < 0  a n d   f → ⋅ L → < 0 d n < 0   a n d   f → ⋅ L → > 0 Fixed d n > 0 d n < 0   a n d   d s > d n tan ⁡ φ d n < 0   a n d   d s < d n tan ⁡ φAs shown in Table1, the contact constraints of the last time step impact the detection of contact constraints of the following time step. For example, if the contact state is open in the last time step, the contact state can be sliding in the following time step only when the equation dn<0 and ds>dntan⁡φ is satisfied. If the contact state is sliding in the last time step, the contact state can be sliding in the following time step only when the equation dn<0 and f→·L→<0 is satisfied. Thus, the type of contact of the last time step must be considered when determining the type of contact in the current time step.


### 3.1. Approximated Step FunctionH1Δ

As shown in Table1, when dn is larger than 0, the type of contact in the current time step will be open regardless of the type of contact of the last time step. Suppose that there exists one step function H1Δ simulating whether embedding occurs: when Δ is smaller than 0, the value of the function is -1 and embedding does not occur; when Δ is larger than 0, the value of the function is 1 and the embedding occurs. The step function H1Δ can be expressed as(5)H1Δ=1,Δ>0-1,Δ<0,where Δ reflects the degree of embedding and is expressed as (1). The step function H1Δ is illustrated in Figure 2.Figure 2
Schematic diagram of embedding simulated by a step function.The specific value ofΔ is not known before the calculation; thus, the specific value of H1Δ cannot be determined. In this case, the concept of an approximated function [25] can be introduced. The following four approximated functions [26, 27] can be observed before approximating H1Δ; they are able to fully approximate the step function, whose value ranges from -1 to 1:(1)
Square root functionG1Δ=Δ/Δ2+ε(2)
Hyperbolic tangent functionG2Δ=tanh⁡Δ/ε(3)
Arctangent functionG3Δ=2/πarctan⁡Δ/ε(4)
Exponential functionG4Δ=1-2/eΔ/ε+1ε is the approximation factor. When ε approaches 0+, the degree of approximation is extremely high. The curves of the three approximated step functions, with ε being 5.0×10-8, are illustrated in Figure 3.Figure 3
Approximated function curve of the step function.In Figure3, the horizontal coordinate Δ reflects the degree of embedding, and the vertical coordinate GΔ reflects the value of the approximated step function. In general, the allowable embedding distance is extremely short; in this case, the magnitude of the horizontal coordinate Δ approaches 10−7, which can meet the computational requirements. The approximation factor ε can be reduced to further improve the accuracy. Although the values of Δ and ε are small, they can be calculated within the existing computational capabilities; thus, the operability and accuracy can be ensured.Figure3 illustrates that the degree of approximation of hyperbolic tangent function G2Δ=tanh⁡Δ/ε is high. When Δ ranges from 0 to 10−7, the slope of the curve of function G2 is high and the value grows rapidly. When Δ reaches 10−7, the curve of function G2 remains fairly flat. From this aspect, the degree of approximation of function G2 is the highest. In contrast, the value of the approximated function grows slightly more slowly when the horizontal coordinate approaches 0, and the degrees of approximation of the other three approximated functions, G1Δ, G3Δ, and G4Δ, are rather low compared to that of G2Δ. The expression of the hyperbolic tangent function G2Δ is relatively simple to express and easy to calculate. Thus, the hyperbolic tangent function G2Δ is used to approximate the step function H1Δ, which can be written as(6)H1Δ=tanh⁡Δε.To ensure that (6) can fully approximate the step function H1Δ, the approximation factor ε in (6) should be sufficiently small and the magnitude should approach 10−8.Thus,H1Δ can be used to simulate whether the current type of contact constraint is open or embedding.


### 3.2. Approximated Step FunctionH2ds/dn

Similarly, embedding is known to occur in the current time step whendn is smaller than 0; then, the type of embedding must be determined. The type of contact of the last time step must be considered. The entire process can be divided into two cases.When the type of contact of the last time step is open or fixed, another step functionH2ds/dn can be used to simulate the fixed embedding and sliding embedding, in which ds is the tangential embedding distance and dn is the normal embedding distance. When calculating tangential embedding distance, we cannot determine the sign of ds, only the absolute value. Suppose that when ds/dn is larger than tan⁡φ or smaller than -tan⁡φ, the value of the function is 1 and sliding embedding occurs; when ds/dn is larger than -tan⁡φ and smaller than tan⁡φ, the value of the function is -1 and fixed embedding occurs, where φ is the friction angle. The function H2ds/dn is illustrated in Figure 4 and it can be expressed as(7)H2dsdn=1,dsdn>tan⁡φordsdn<-tan⁡φ-1,-tan⁡φ<dsdn<tan⁡φ.With reference to (6), H2ds/dn can be approximated using the following equation:(8)H2dsdn=tanh⁡ds/dn-tan⁡φε-tanhds/dn+tan⁡φε+1.Figure 4
Schematic diagram of the embedding mode simulated by a step function (open or fixed in the last calculation step).Thus,H2ds/dn can be used to simulate whether the current state of the contact constraint is fixed embedding or sliding embedding when the type of contact of the last time step is open or fixed.


### 3.3. Approximated Step FunctionH3f→·L→

When the type of contact of the last time step is sliding, another step functionH3f→·L→ can be used to simulate the fixed embedding and sliding embedding, in which f→ is the friction vector and L→ is the shear displacement vector. The friction vector f→ is the product of the tangential stiffness K and the shear displacement vector L′⃗ of the last time step, and the shear displacement vector L→ of the current time step is the vector of the tangential embedding distance ds in Section 2. Suppose that when f→·L→ is larger than 0, the value of the function is 1 and sliding embedding occurs, whereas when ds/dn is smaller than 0, the value of the function is -1 and fixed embedding occurs. The function H3f→·L→ is illustrated in Figure 5 and it can be expressed as(9)H3f→·L→=1,f→·L→>0-1,f→·L→<0.With reference to (6), H3f→·L→ can be approximated using the following equation:(10)H3f→·L→=tanh⁡f→·L→ε.Figure 5
Schematic diagram of the embedding mode simulated by a step function (sliding in the last calculation step).Thus,H3f→·L→ can be used to simulate whether the current state of the contact constraint is fixed embedding or sliding embedding when the type of contact of the last time step is sliding.The variablesΔ, ds, dn, and L→ are functions containing only displacement variables D as an unknown; thus, (6), (8), and (10) are also functions containing only displacement variables D as an unknown; the three equations above can reflect contact constraints (open, sliding, and fixed). From this perspective, the relationship between displacement variables and contact constraints between blocks C=gD can be generalized using H1Δ, H2ds/dn, and H3f→·L→. Next, the concept of Lagrange interpolation is introduced to deduce the expression of the contact potential energy function Πc.


## 3.1. Approximated Step FunctionH1Δ

As shown in Table1, when dn is larger than 0, the type of contact in the current time step will be open regardless of the type of contact of the last time step. Suppose that there exists one step function H1Δ simulating whether embedding occurs: when Δ is smaller than 0, the value of the function is -1 and embedding does not occur; when Δ is larger than 0, the value of the function is 1 and the embedding occurs. The step function H1Δ can be expressed as(5)H1Δ=1,Δ>0-1,Δ<0,where Δ reflects the degree of embedding and is expressed as (1). The step function H1Δ is illustrated in Figure 2.Figure 2
Schematic diagram of embedding simulated by a step function.The specific value ofΔ is not known before the calculation; thus, the specific value of H1Δ cannot be determined. In this case, the concept of an approximated function [25] can be introduced. The following four approximated functions [26, 27] can be observed before approximating H1Δ; they are able to fully approximate the step function, whose value ranges from -1 to 1:(1)
Square root functionG1Δ=Δ/Δ2+ε(2)
Hyperbolic tangent functionG2Δ=tanh⁡Δ/ε(3)
Arctangent functionG3Δ=2/πarctan⁡Δ/ε(4)
Exponential functionG4Δ=1-2/eΔ/ε+1ε is the approximation factor. When ε approaches 0+, the degree of approximation is extremely high. The curves of the three approximated step functions, with ε being 5.0×10-8, are illustrated in Figure 3.Figure 3
Approximated function curve of the step function.In Figure3, the horizontal coordinate Δ reflects the degree of embedding, and the vertical coordinate GΔ reflects the value of the approximated step function. In general, the allowable embedding distance is extremely short; in this case, the magnitude of the horizontal coordinate Δ approaches 10−7, which can meet the computational requirements. The approximation factor ε can be reduced to further improve the accuracy. Although the values of Δ and ε are small, they can be calculated within the existing computational capabilities; thus, the operability and accuracy can be ensured.Figure3 illustrates that the degree of approximation of hyperbolic tangent function G2Δ=tanh⁡Δ/ε is high. When Δ ranges from 0 to 10−7, the slope of the curve of function G2 is high and the value grows rapidly. When Δ reaches 10−7, the curve of function G2 remains fairly flat. From this aspect, the degree of approximation of function G2 is the highest. In contrast, the value of the approximated function grows slightly more slowly when the horizontal coordinate approaches 0, and the degrees of approximation of the other three approximated functions, G1Δ, G3Δ, and G4Δ, are rather low compared to that of G2Δ. The expression of the hyperbolic tangent function G2Δ is relatively simple to express and easy to calculate. Thus, the hyperbolic tangent function G2Δ is used to approximate the step function H1Δ, which can be written as(6)H1Δ=tanh⁡Δε.To ensure that (6) can fully approximate the step function H1Δ, the approximation factor ε in (6) should be sufficiently small and the magnitude should approach 10−8.Thus,H1Δ can be used to simulate whether the current type of contact constraint is open or embedding.


## 3.2. Approximated Step FunctionH2ds/dn

Similarly, embedding is known to occur in the current time step whendn is smaller than 0; then, the type of embedding must be determined. The type of contact of the last time step must be considered. The entire process can be divided into two cases.When the type of contact of the last time step is open or fixed, another step functionH2ds/dn can be used to simulate the fixed embedding and sliding embedding, in which ds is the tangential embedding distance and dn is the normal embedding distance. When calculating tangential embedding distance, we cannot determine the sign of ds, only the absolute value. Suppose that when ds/dn is larger than tan⁡φ or smaller than -tan⁡φ, the value of the function is 1 and sliding embedding occurs; when ds/dn is larger than -tan⁡φ and smaller than tan⁡φ, the value of the function is -1 and fixed embedding occurs, where φ is the friction angle. The function H2ds/dn is illustrated in Figure 4 and it can be expressed as(7)H2dsdn=1,dsdn>tan⁡φordsdn<-tan⁡φ-1,-tan⁡φ<dsdn<tan⁡φ.With reference to (6), H2ds/dn can be approximated using the following equation:(8)H2dsdn=tanh⁡ds/dn-tan⁡φε-tanhds/dn+tan⁡φε+1.Figure 4
Schematic diagram of the embedding mode simulated by a step function (open or fixed in the last calculation step).Thus,H2ds/dn can be used to simulate whether the current state of the contact constraint is fixed embedding or sliding embedding when the type of contact of the last time step is open or fixed.


## 3.3. Approximated Step FunctionH3f→·L→

When the type of contact of the last time step is sliding, another step functionH3f→·L→ can be used to simulate the fixed embedding and sliding embedding, in which f→ is the friction vector and L→ is the shear displacement vector. The friction vector f→ is the product of the tangential stiffness K and the shear displacement vector L′⃗ of the last time step, and the shear displacement vector L→ of the current time step is the vector of the tangential embedding distance ds in Section 2. Suppose that when f→·L→ is larger than 0, the value of the function is 1 and sliding embedding occurs, whereas when ds/dn is smaller than 0, the value of the function is -1 and fixed embedding occurs. The function H3f→·L→ is illustrated in Figure 5 and it can be expressed as(9)H3f→·L→=1,f→·L→>0-1,f→·L→<0.With reference to (6), H3f→·L→ can be approximated using the following equation:(10)H3f→·L→=tanh⁡f→·L→ε.Figure 5
Schematic diagram of the embedding mode simulated by a step function (sliding in the last calculation step).Thus,H3f→·L→ can be used to simulate whether the current state of the contact constraint is fixed embedding or sliding embedding when the type of contact of the last time step is sliding.The variablesΔ, ds, dn, and L→ are functions containing only displacement variables D as an unknown; thus, (6), (8), and (10) are also functions containing only displacement variables D as an unknown; the three equations above can reflect contact constraints (open, sliding, and fixed). From this perspective, the relationship between displacement variables and contact constraints between blocks C=gD can be generalized using H1Δ, H2ds/dn, and H3f→·L→. Next, the concept of Lagrange interpolation is introduced to deduce the expression of the contact potential energy function Πc.


## 4. Contact Potential Energy Function Based on Lagrange Interpolation

Approximated step functionsH1Δ, H2ds/dn, and H3f→·L→ can yield a certain value (-1 to 1), which can be directly calculated to obtain the contact constraints without open-close iteration. The rules of contact detection with an approximated step function are shown in Table 2.Table 2
Discriminant table of contact constraints with an approximated step function.
Time step →i-1 Time step →i Open Sliding Fixed Open or fixed H 1 Δ = - 1 H 1 Δ = 1   a n d   H 2 d s d n = 1 H 1 Δ = 1   a n d   H 2 d s d n = - 1 Sliding H 1 Δ = - 1 H 1 Δ = 1   a n d   H 3 f → ⋅ L → = 1 H 1 Δ = 1    a n d   H 3 f → ⋅ L → = - 1The potential energy function in response to different types of contact constraints is not the same but has a one-to-one correspondence. When fixed embedding occurs, there exists a fixed contact potential energy functionΠfix; when sliding embedding occurs, there exists a sliding contact potential energy function Πslide. Based on Lagrange interpolation and Table 2, the contact potential energy function can be expressed as detailed below.When the type of contact of the last time step is open or fixed, it can be written as(11)Πc-i=H1Δ--11--1×H2ds/dn--11--1×Πslide-i+H1Δ--11--1×H2ds/dn-1-1-1×Πfix-i+H1Δ-1-1-1×0i=1,…,k.When the type of contact of the last time step is sliding, it can be written as(12)Πc-i=H1Δ--11--1×H3f→·L→--11--1×Πslide-i+H1Δ--11--1×H3f→·L→-1-1-1×Πfix-i+H1Δ-1-1-1×0i=1,…,k.Equations (11) and (12) represent the current contact potential energy function with different types of contact constraints of the last time step, where k is the number of contact pairs. The fixed contact potential energy Πfix and sliding contact potential energy Πfix are based on the conventional DDA method [1] and contain only displacement variables D as an unknown.Suppose that the number of contact pairs isk. Then, (13) can be deduced by summing up all the contact potential energy values.(13)Πcsum=∑i=1kΠc-i.Equation (13) represents the overall contact potential energy Πcsum of the block system. The potential energy inside blocks Πbsum is generated with stress, strain, and loading. For details regarding the block potential energy Πbsum, refer to the conventional DDA method [1]. The overall potential energy can be deduced by summing up the contact potential energy Πcsum and the block potential energy Πbsum as(14)Πsum=Πcsum+Πbsum.The contact constraints have already been represented using the displacement of blocks, and (14) contains only the displacement variables D as an unknown(15)Πsum=ΠsumD.When taking the derivative of (15), there is no need to consider the contact constraints and the implicit open-close iteration.Next, the variable metric method is used in (15) to analyze the extremum of the potential energy function of the block system.


## 5. Minimization Problem of the Overall Potential Energy Solved by the Variable Metric Method

The variable metric method is an efficient algorithm for solving the extreme problem of an unconstrained and multidimensional function and is widely applied in geotechnical engineering [28]. This method does not require the calculation of the second derivative matrix and its inverse operation. Its convergence rate is considerably faster than that of gradient method’s and exhibits second-order global convergence, including high-dimensional problems. Thus, the variable metric method is chosen to determine the overall potential energy of a block system with the improved DDA method based on explicitly solving contact constraints.


### 5.1. Iteration Scheme of the Variable Metric Method

The variable metric method is developed based on Newton’s method; the basic iteration scheme constructed can be written as follows:(16)Dk+1=Dk-αkAk∇ΠDk,where Dk+1 and Dk are the displacement variables of blocks at time steps k+1 and k, respectively, Π is the overall potential energy, ∇Π is the gradient of the overall potential energy with respect to the displacement variables of blocks D, α(k) is the optimal step length of time step k, and A(k) is a symmetric matrix with size N×N (where N is the number of freedom of the overall displacement variables). Next, the iteration factors mentioned above are discussed in detail.


### 5.2. Determination of Iteration Factors

The three main iteration factors are the gradient of the overall potential energy∇Π, the optimal step length factor α(k), and the symmetric matrix A(k).Currently, the ahead difference quotient is widely used to approximate the gradient, with∇Π¯Dk replacing ∇ΠDk; the expression can be written as(17)∇Π¯Dk=ΠDk+hke-ΠDkhk,where hk controls the accuracy of the calculation and e is the unit vector in n-dimensional Euclidean space.α ( k ) is the optimal step length factor of time step k. The basic process of realizing this method is determining the value of the optimal step length αk with a certain descent direction, which could lead to a sufficient decrease in the value of the objective function. The most widely used methods are the Wolfe method and Goldstein method [29].The variable metric matrixA(k) can fully approximate the inverse matrix of the Hesse matrix ∇2ΠDk-1 when reaching the extremum. The variable metric matrix of the adjacent time step has the following relationship:(18)Ak+1=Ak+ΔAk,where ΔAk is the calibrated matrix and can be obtained using the DFP method, written as follows:(19)ΔAk=ΔDkΔDkTΔDkTΔgk-AkΔgkΔgkTAkΔgkTAkΔgk,where Dk is the displacement variables of blocks at time step k, ΔDk is the displacement difference ΔDk=Dk+1-Dk, gk is the gradient of the overall potential energy of the displacement of blocks at time step k, and Δgk is the gradient difference Δgk=gk+1-gk.


### 5.3. Process of Iteration

The optimal solution of displacement variables is obtained with the following steps:(1)
Determine the overall potential energy function∇ΠD deduced in Sections 3 and 4.(2)
Specify the initial iteration pointD0 and termination factor ε′.(3)
Set time stepk being 0 and A0 being identity matrix I.(4)
Calculate the value of∇ΠDk referred to in (17) and set ∇ΠDk·Ak being Sk.(5)
Set the value ofαk in order to minimize the value of ΠDk-αkSk which can be referred to in [29].(6)
Check if the termination condition∇ΠDk≤ε′ is verified. If so, stop iterating and set Dk+1=Dk-αkSk being the optimal solution of displacement variables. If not, conduct the following step.(7)
Calculate the value ofΔAk and Ak+1 referred to in (18) and (19). Set time step k being k+1 and turn to step 4.After we obtain the displacement of blocksD, the contact constraints can be obtained by substituting D into C=gD. The relationship between displacement variables and contact constraints between blocks C=gD has been generalized using H1Δ, H2ds/dn, and H3f→·L→ in Section 3 and there is no need to conduct an implicit open-close iteration.


## 5.1. Iteration Scheme of the Variable Metric Method

The variable metric method is developed based on Newton’s method; the basic iteration scheme constructed can be written as follows:(16)Dk+1=Dk-αkAk∇ΠDk,where Dk+1 and Dk are the displacement variables of blocks at time steps k+1 and k, respectively, Π is the overall potential energy, ∇Π is the gradient of the overall potential energy with respect to the displacement variables of blocks D, α(k) is the optimal step length of time step k, and A(k) is a symmetric matrix with size N×N (where N is the number of freedom of the overall displacement variables). Next, the iteration factors mentioned above are discussed in detail.


## 5.2. Determination of Iteration Factors

The three main iteration factors are the gradient of the overall potential energy∇Π, the optimal step length factor α(k), and the symmetric matrix A(k).Currently, the ahead difference quotient is widely used to approximate the gradient, with∇Π¯Dk replacing ∇ΠDk; the expression can be written as(17)∇Π¯Dk=ΠDk+hke-ΠDkhk,where hk controls the accuracy of the calculation and e is the unit vector in n-dimensional Euclidean space.α ( k ) is the optimal step length factor of time step k. The basic process of realizing this method is determining the value of the optimal step length αk with a certain descent direction, which could lead to a sufficient decrease in the value of the objective function. The most widely used methods are the Wolfe method and Goldstein method [29].The variable metric matrixA(k) can fully approximate the inverse matrix of the Hesse matrix ∇2ΠDk-1 when reaching the extremum. The variable metric matrix of the adjacent time step has the following relationship:(18)Ak+1=Ak+ΔAk,where ΔAk is the calibrated matrix and can be obtained using the DFP method, written as follows:(19)ΔAk=ΔDkΔDkTΔDkTΔgk-AkΔgkΔgkTAkΔgkTAkΔgk,where Dk is the displacement variables of blocks at time step k, ΔDk is the displacement difference ΔDk=Dk+1-Dk, gk is the gradient of the overall potential energy of the displacement of blocks at time step k, and Δgk is the gradient difference Δgk=gk+1-gk.


## 5.3. Process of Iteration

The optimal solution of displacement variables is obtained with the following steps:(1)
Determine the overall potential energy function∇ΠD deduced in Sections 3 and 4.(2)
Specify the initial iteration pointD0 and termination factor ε′.(3)
Set time stepk being 0 and A0 being identity matrix I.(4)
Calculate the value of∇ΠDk referred to in (17) and set ∇ΠDk·Ak being Sk.(5)
Set the value ofαk in order to minimize the value of ΠDk-αkSk which can be referred to in [29].(6)
Check if the termination condition∇ΠDk≤ε′ is verified. If so, stop iterating and set Dk+1=Dk-αkSk being the optimal solution of displacement variables. If not, conduct the following step.(7)
Calculate the value ofΔAk and Ak+1 referred to in (18) and (19). Set time step k being k+1 and turn to step 4.After we obtain the displacement of blocksD, the contact constraints can be obtained by substituting D into C=gD. The relationship between displacement variables and contact constraints between blocks C=gD has been generalized using H1Δ, H2ds/dn, and H3f→·L→ in Section 3 and there is no need to conduct an implicit open-close iteration.


## 6. Implementation Procedure

The difference in contact detection between the conventional DDA method and improved DDA method is illustrated in Figure6.Figure 6
Comparison of contact detection between the two DDA methods.The contact constraints are simulated by approximated step function and Lagrange interpolation. There is no need to conduct implicit open-close iteration in the improved DDA method. The calculation process is illustrated in Figure7.Figure 7
Calculation process of the improved DDA method.


## 7. Validation of the Improved DDA Method

Based on the improved DDA method proposed in this paper, the authors have written the corresponding calculation program code with Fortran and applied it to the following two numerical models. A model of a sliding block is used to test the accuracy of the improved DDA method, and a model of a tunnel is used to test the computational efficiency and robustness.


### 7.1. Model of a Sliding Block

#### 7.1.1. Numerical Model

To validate the accuracy and iterative stability of the proposed method, a simple but representative model of a sliding block is used to compare the numerical results and analytical results.The model of a sliding block is illustrated in Figure8; it consists of a free sliding block and a base of variable slope [30]. The bottom and lateral sides of the base are fixed, and displacement cannot occur. The block can slide along the slope without any constraint. The dip angle θ1 of the upper part of the slope is 50°, and the length L1 is 8 m. The dip angle θ2 of the lower part of the slope is 20°, and the length L2 is 18 m. Young’s modulus E of the block and base is 12 MPa, the Poisson ratio is 0.32, and the value of the contact stiffness is 60 times that of Young’s modulus.Figure 8
Model of a sliding block on a slope with a variable gradient.To minimize the energy loss at the turning point, the size of the block should be sufficiently small for the transition to be smooth. The block weighs 500 kg, with dimensions of 0.6 × 1.2 m.Because the model of a sliding block is rather simple, it is easy to analyze the analytical solution of the block over time.


#### 7.1.2. Calculation Results

Suppose that in this example the cohesive force on the sliding surface is 0. The motion of the sliding block is analyzed using a friction angleφ of 30°. Theoretically, the friction angle φ is smaller than the dip angle of the upper part of slope θ1 but larger than the dip angle of the lower part θ2; thus, the sliding block would first accelerate to slide and then slow down until stopping. The displacement parallel to the sliding surface is calculated using the theoretical method, the conventional DDA method, and the improved DDA method. The absolute error of the theoretical results s′ and the two types of DDA results s are analyzed. r=s-s′/s is taken as the relative error. The results are shown in Figures 9–11.Figure 9
Absolute error-time curve between two DDA solutions and analytical solution.Figure 10
Relative error-time curve between two DDA solutions and analytical solution (overview).Figure 11
Relative error-time curve between two DDA solutions and analytical solution (partially zoomed-in view).As shown in Figure9, the results of the variation law of the two curves of DDA are nearly the same, but there are subtle differences in values. In the early stage of movement, the block has not reached a steady state of motion, and the contact constraints between the block and base vary from no contact to slight embedding, which leads to irregular displacement and a certain absolute error. As time progresses, the error accumulates and the absolute error grows, but the state of motion remains steady. Comparing the results of the two DDA methods illustrates that the absolute error of the results of the improved DDA method is slightly larger than that of the conventional DDA method.Because the absolute error is relatively small compared with the actual displacement and the absolute errors of the two DDA methods in Figure9 are similar, the curves of the relative error of the two DDA methods in Figure 10 largely coincide. The relative error in the early stage reaches a maximum, with a value of 3.1%. As the motion remains steady, the relative error grows slightly and ultimately approaches 0.The relative error in Figure10 spans over a long period of time. In the early stage of motion, the relative error is relatively large and the curves are relatively steep. At a later stage of motion, the relative error approaches 0, and the curves are relatively flat. The curves of the relative error of the two DDA methods largely coincide, and there exists homogeneous slope variation in the two curves only for the time span of 0.1–1.0 s. The curves of the relative error of the two DDA methods for this time span are shown in Figure 11; the variations in the slope are notable. The relative error of the improved DDA method is slightly larger than that of the conventional DDA method.As illustrated in Figures9, 10, and 11, in the early stage of motion, both the absolute error and relative error of the improved DDA method are slightly larger than those of the conventional DDA method. The contact constraints can be precisely determined with open-close iteration, but the convergence of the open-close iteration is poor. The improved DDA method can simulate the contact constraints with the step function approximated using the hyperbolic tangent function; the calculation accuracy is related to the degree of embedding Δ and the approximation factor ε. The calculating accuracies of the two DDA methods illustrate that the error ranges of the two are highly similar. The difference value of the two DDA methods regarding the absolute error is largest at the beginning of motion, reaching 0.02 mm, and the difference value of the two DDA methods regarding relative error is no more than 0.2% when it reaches the maximum. At a later stage of motion, the error curves of the two DDA methods largely coincide. The improved DDA method can control the calculation accuracy with the approximation factor ε and the termination condition ∇ΠDk≤ε′. Thus, the improved DDA method based on explicitly solving contact constraints has a high calculation accuracy and can satisfy the engineering calculation requirements.Contact force is generated when the block is sliding on the slope, and it has a direct impact on the displacement and calculating accuracy of the block [31]. When the friction angle φ is 54°, the block remains at the initial position, and the normal and tangential embedding distances of the two vertexes at the bottom can be calculated to validate the calculation accuracy of the contact force. The theoretical embedding distance [32] could be obtained using (20), which is written as follows:(20)DNL=Wcos⁡θ0.5+h/2btan⁡θK,DNU=Wcos⁡θ0.5-h/2btan⁡θK,DS=Wsin⁡θK,where DNL and DNU are the normal embedding distances of the lower contact point PL and upper contact point PU, respectively, DS is the sum of all the tangential embedding distances, W is the weight of the block, h and b are the height and width of the block, respectively, and θ is the dip angle of the slope.Comparing the analytical resultd and DDA result d′, r=d-d′/d is taken as the relative error; the results are shown in Table 3.Table 3
Comparison of the embedding distances for block stabilized using the three calculation solutions.
Contact point Analytical solution Conventional DDA method solution Improved DDA method solution Embedding distance (10−4 m) Embedding distance (10−4 m) Relative error (%) Embedding distance (10−4 m) Relative error (%) DNL 2.0944 2.0759 0.8833 2.0684 1.1214 DNU 0.5304 0.5254 0.9426 0.5236 1.2255 DSL — 1.8186 — 1.8097 — DSU — 1.3094 — 1.3086 — DS 3.1280 3.1016 0.8439 3.0894 0.9463As shown in Table3, the difference between the DDA results and analytical results is small and the relative error is controlled within 1.5%, which illustrates the proposed method’s high accuracy. Thus, it can be assumed that the calculating accuracy of the contact force is high and can satisfy the engineering calculation requirements.


### 7.2. Model of a Tunnel

#### 7.2.1. Numerical Model

The model of an underground tunnel with fully developed joints is shown in Figure12; the rock mass is cut into blocks by the joint planes. The blocks have the following material properties: bulk density γ=27KN/m3, Young’s modulus E=1.50 MPa, Poisson’s ratio μ=0.23, friction angle φ=5°, cohesive force c=0, and contact stiffness P=50E. All blocks are subjected to the same initial stresses of σx=-0.05 MPa, σy=-1.05 MPa, and τxy=0, and there is no body force.Figure 12
Model of the tunnel before deformation.


#### 7.2.2. Calculation Results

The time step is taken asΔ=0.01 s, and the deformation of the tunnel after 100 time steps is illustrated in Figure 13. The surrounding rock is in a state of self-stabilization, and the deformation and displacement of the block in the four corners are large compared with those in other locations. The outer contour of the tunnel is fixed, and no displacement occurs.Figure 13
Model of the tunnel after deformation.The conventional DDA method and improved DDA method differ regarding the calculation process. The difference mainly concerns the contact detection and solving of the extremum of the potential energy, not the pretreatment, postprocessing, and analysis of the potential energy inside the blocks. In the conventional DDA method, implicit open-close iteration and solving the large sparse linear equations are required at every time step. Implicit open-close iteration can be viewed as the process of determining the contact constraint, whereas solving the large sparse linear equations can be seen as the process of analyzing the extremum of the potential energy. More time is consumed conducting implicit open-close iteration compared with solving large sparse linear equations. In contrast, in the improved DDA method, there is no need to conduct implicit open-close iteration, and the contact constraints between blocks have already been expressed using the approximated step function and Lagrange interpolation. The expression of the overall potential energy function of the improved DDA method is more complex than that of the conventional DDA method, and the time is mainly consumed while determining the extremum of the potential energy.Because the overall potential energy function and the process of solving equations differ between the two DDA methods, the iteration step cannot be used to evaluate the efficiency of the two methods. In this paper, the time consumed during calculation is used as the comparison criterion. Because the pretreatment, postprocessing, and analysis of the potential energy inside blocks of the two DDA methods are the same and their difference with respect to time consumed is mainly caused by contact detection and solving the extremum of the potential energy, it is reasonable to take the time consumed as the criterion. All analyses are performed on the same computer with the following system configuration: processor speed = 3.40 GHz and RAM = 8.0 GB.Four calculation conditions, with contact stiffness of 80E, 100E, 120E, and 150E, are given in Table 4. The conventional DDA method and improved DDA method are used to analyze the stability of the surrounding rock of the tunnel after 100 time steps, and the time step is taken as Δ=0.01 s. The time consumed is recorded, and the results are shown in Table 4.Table 4
Comparison of iteration time with different contact stiffness.
Contact stiffness Iteration time (s) Conventional DDA method Improved DDA method 5 0 E 793.2 592.6 80 E 674.9 602.4 100 E 812.7 613.7 120 E Terminated 589.5 150 E 931.9 598.3As shown in Table4, the time consumed and the variation trends of the two DDA methods are different. The improved DDA method is more stable, and the time consumed for different contact stiffness remains at approximately 600 s; the contact stiffness has only a slight impact on the time consumed. Compared with the improved DDA method, the time consumed using the conventional DDA method is generally over 600 s, the variation trends are not distinct, and the embedding distance varies with changing contact stiffness. Because the contact constraint is uncertain in the conventional DDA method, considerable time is spent performing the implicit open-close iteration. If the open-close iteration is successful, it needs to be conducted only a few times to obtain the precise contact constraints; thus, the amount of time consumed is not overly burdensome. If the precise contact constraint cannot be obtained after conducting the implicit open-close iteration more than six times, the time step will be reduced to 1/3 of the original. In this case, the time step will be smaller and considerable time will be spent in the previous process; thus, considerably more time is consumed by the conventional DDA method than by the improved DDA method.When the value of the contact stiffness is 120E, the calculation process with the conventional DDA method cannot be completed. The sign “terminated” means that, in the process of implicit open-close iteration, the precise contact constraints could not be found and the time step was reduced several times, finally falling into an infinite loop and being forced to terminate [33]. The robustness of the conventional DDA method is not strong, and in some cases, the implicit open-close iteration cannot be conducted successfully, thus preventing the final correct solution from being obtained. Because implicit open-close iteration is not needed in the improved DDA method, the above problems do not occur. Thus, the improved DDA method has a more stable and more robust calculation convergence.Based on the two aspects above, the improved DDA method exhibits a faster calculation speed than the conventional DDA method; moreover, it has a more robust calculation convergence.


## 7.1. Model of a Sliding Block

### 7.1.1. Numerical Model

To validate the accuracy and iterative stability of the proposed method, a simple but representative model of a sliding block is used to compare the numerical results and analytical results.The model of a sliding block is illustrated in Figure8; it consists of a free sliding block and a base of variable slope [30]. The bottom and lateral sides of the base are fixed, and displacement cannot occur. The block can slide along the slope without any constraint. The dip angle θ1 of the upper part of the slope is 50°, and the length L1 is 8 m. The dip angle θ2 of the lower part of the slope is 20°, and the length L2 is 18 m. Young’s modulus E of the block and base is 12 MPa, the Poisson ratio is 0.32, and the value of the contact stiffness is 60 times that of Young’s modulus.Figure 8
Model of a sliding block on a slope with a variable gradient.To minimize the energy loss at the turning point, the size of the block should be sufficiently small for the transition to be smooth. The block weighs 500 kg, with dimensions of 0.6 × 1.2 m.Because the model of a sliding block is rather simple, it is easy to analyze the analytical solution of the block over time.


### 7.1.2. Calculation Results

Suppose that in this example the cohesive force on the sliding surface is 0. The motion of the sliding block is analyzed using a friction angleφ of 30°. Theoretically, the friction angle φ is smaller than the dip angle of the upper part of slope θ1 but larger than the dip angle of the lower part θ2; thus, the sliding block would first accelerate to slide and then slow down until stopping. The displacement parallel to the sliding surface is calculated using the theoretical method, the conventional DDA method, and the improved DDA method. The absolute error of the theoretical results s′ and the two types of DDA results s are analyzed. r=s-s′/s is taken as the relative error. The results are shown in Figures 9–11.Figure 9
Absolute error-time curve between two DDA solutions and analytical solution.Figure 10
Relative error-time curve between two DDA solutions and analytical solution (overview).Figure 11
Relative error-time curve between two DDA solutions and analytical solution (partially zoomed-in view).As shown in Figure9, the results of the variation law of the two curves of DDA are nearly the same, but there are subtle differences in values. In the early stage of movement, the block has not reached a steady state of motion, and the contact constraints between the block and base vary from no contact to slight embedding, which leads to irregular displacement and a certain absolute error. As time progresses, the error accumulates and the absolute error grows, but the state of motion remains steady. Comparing the results of the two DDA methods illustrates that the absolute error of the results of the improved DDA method is slightly larger than that of the conventional DDA method.Because the absolute error is relatively small compared with the actual displacement and the absolute errors of the two DDA methods in Figure9 are similar, the curves of the relative error of the two DDA methods in Figure 10 largely coincide. The relative error in the early stage reaches a maximum, with a value of 3.1%. As the motion remains steady, the relative error grows slightly and ultimately approaches 0.The relative error in Figure10 spans over a long period of time. In the early stage of motion, the relative error is relatively large and the curves are relatively steep. At a later stage of motion, the relative error approaches 0, and the curves are relatively flat. The curves of the relative error of the two DDA methods largely coincide, and there exists homogeneous slope variation in the two curves only for the time span of 0.1–1.0 s. The curves of the relative error of the two DDA methods for this time span are shown in Figure 11; the variations in the slope are notable. The relative error of the improved DDA method is slightly larger than that of the conventional DDA method.As illustrated in Figures9, 10, and 11, in the early stage of motion, both the absolute error and relative error of the improved DDA method are slightly larger than those of the conventional DDA method. The contact constraints can be precisely determined with open-close iteration, but the convergence of the open-close iteration is poor. The improved DDA method can simulate the contact constraints with the step function approximated using the hyperbolic tangent function; the calculation accuracy is related to the degree of embedding Δ and the approximation factor ε. The calculating accuracies of the two DDA methods illustrate that the error ranges of the two are highly similar. The difference value of the two DDA methods regarding the absolute error is largest at the beginning of motion, reaching 0.02 mm, and the difference value of the two DDA methods regarding relative error is no more than 0.2% when it reaches the maximum. At a later stage of motion, the error curves of the two DDA methods largely coincide. The improved DDA method can control the calculation accuracy with the approximation factor ε and the termination condition ∇ΠDk≤ε′. Thus, the improved DDA method based on explicitly solving contact constraints has a high calculation accuracy and can satisfy the engineering calculation requirements.Contact force is generated when the block is sliding on the slope, and it has a direct impact on the displacement and calculating accuracy of the block [31]. When the friction angle φ is 54°, the block remains at the initial position, and the normal and tangential embedding distances of the two vertexes at the bottom can be calculated to validate the calculation accuracy of the contact force. The theoretical embedding distance [32] could be obtained using (20), which is written as follows:(20)DNL=Wcos⁡θ0.5+h/2btan⁡θK,DNU=Wcos⁡θ0.5-h/2btan⁡θK,DS=Wsin⁡θK,where DNL and DNU are the normal embedding distances of the lower contact point PL and upper contact point PU, respectively, DS is the sum of all the tangential embedding distances, W is the weight of the block, h and b are the height and width of the block, respectively, and θ is the dip angle of the slope.Comparing the analytical resultd and DDA result d′, r=d-d′/d is taken as the relative error; the results are shown in Table 3.Table 3
Comparison of the embedding distances for block stabilized using the three calculation solutions.
Contact point Analytical solution Conventional DDA method solution Improved DDA method solution Embedding distance (10−4 m) Embedding distance (10−4 m) Relative error (%) Embedding distance (10−4 m) Relative error (%) DNL 2.0944 2.0759 0.8833 2.0684 1.1214 DNU 0.5304 0.5254 0.9426 0.5236 1.2255 DSL — 1.8186 — 1.8097 — DSU — 1.3094 — 1.3086 — DS 3.1280 3.1016 0.8439 3.0894 0.9463As shown in Table3, the difference between the DDA results and analytical results is small and the relative error is controlled within 1.5%, which illustrates the proposed method’s high accuracy. Thus, it can be assumed that the calculating accuracy of the contact force is high and can satisfy the engineering calculation requirements.


## 7.1.1. Numerical Model

To validate the accuracy and iterative stability of the proposed method, a simple but representative model of a sliding block is used to compare the numerical results and analytical results.The model of a sliding block is illustrated in Figure8; it consists of a free sliding block and a base of variable slope [30]. The bottom and lateral sides of the base are fixed, and displacement cannot occur. The block can slide along the slope without any constraint. The dip angle θ1 of the upper part of the slope is 50°, and the length L1 is 8 m. The dip angle θ2 of the lower part of the slope is 20°, and the length L2 is 18 m. Young’s modulus E of the block and base is 12 MPa, the Poisson ratio is 0.32, and the value of the contact stiffness is 60 times that of Young’s modulus.Figure 8
Model of a sliding block on a slope with a variable gradient.To minimize the energy loss at the turning point, the size of the block should be sufficiently small for the transition to be smooth. The block weighs 500 kg, with dimensions of 0.6 × 1.2 m.Because the model of a sliding block is rather simple, it is easy to analyze the analytical solution of the block over time.


## 7.1.2. Calculation Results

Suppose that in this example the cohesive force on the sliding surface is 0. The motion of the sliding block is analyzed using a friction angleφ of 30°. Theoretically, the friction angle φ is smaller than the dip angle of the upper part of slope θ1 but larger than the dip angle of the lower part θ2; thus, the sliding block would first accelerate to slide and then slow down until stopping. The displacement parallel to the sliding surface is calculated using the theoretical method, the conventional DDA method, and the improved DDA method. The absolute error of the theoretical results s′ and the two types of DDA results s are analyzed. r=s-s′/s is taken as the relative error. The results are shown in Figures 9–11.Figure 9
Absolute error-time curve between two DDA solutions and analytical solution.Figure 10
Relative error-time curve between two DDA solutions and analytical solution (overview).Figure 11
Relative error-time curve between two DDA solutions and analytical solution (partially zoomed-in view).As shown in Figure9, the results of the variation law of the two curves of DDA are nearly the same, but there are subtle differences in values. In the early stage of movement, the block has not reached a steady state of motion, and the contact constraints between the block and base vary from no contact to slight embedding, which leads to irregular displacement and a certain absolute error. As time progresses, the error accumulates and the absolute error grows, but the state of motion remains steady. Comparing the results of the two DDA methods illustrates that the absolute error of the results of the improved DDA method is slightly larger than that of the conventional DDA method.Because the absolute error is relatively small compared with the actual displacement and the absolute errors of the two DDA methods in Figure9 are similar, the curves of the relative error of the two DDA methods in Figure 10 largely coincide. The relative error in the early stage reaches a maximum, with a value of 3.1%. As the motion remains steady, the relative error grows slightly and ultimately approaches 0.The relative error in Figure10 spans over a long period of time. In the early stage of motion, the relative error is relatively large and the curves are relatively steep. At a later stage of motion, the relative error approaches 0, and the curves are relatively flat. The curves of the relative error of the two DDA methods largely coincide, and there exists homogeneous slope variation in the two curves only for the time span of 0.1–1.0 s. The curves of the relative error of the two DDA methods for this time span are shown in Figure 11; the variations in the slope are notable. The relative error of the improved DDA method is slightly larger than that of the conventional DDA method.As illustrated in Figures9, 10, and 11, in the early stage of motion, both the absolute error and relative error of the improved DDA method are slightly larger than those of the conventional DDA method. The contact constraints can be precisely determined with open-close iteration, but the convergence of the open-close iteration is poor. The improved DDA method can simulate the contact constraints with the step function approximated using the hyperbolic tangent function; the calculation accuracy is related to the degree of embedding Δ and the approximation factor ε. The calculating accuracies of the two DDA methods illustrate that the error ranges of the two are highly similar. The difference value of the two DDA methods regarding the absolute error is largest at the beginning of motion, reaching 0.02 mm, and the difference value of the two DDA methods regarding relative error is no more than 0.2% when it reaches the maximum. At a later stage of motion, the error curves of the two DDA methods largely coincide. The improved DDA method can control the calculation accuracy with the approximation factor ε and the termination condition ∇ΠDk≤ε′. Thus, the improved DDA method based on explicitly solving contact constraints has a high calculation accuracy and can satisfy the engineering calculation requirements.Contact force is generated when the block is sliding on the slope, and it has a direct impact on the displacement and calculating accuracy of the block [31]. When the friction angle φ is 54°, the block remains at the initial position, and the normal and tangential embedding distances of the two vertexes at the bottom can be calculated to validate the calculation accuracy of the contact force. The theoretical embedding distance [32] could be obtained using (20), which is written as follows:(20)DNL=Wcos⁡θ0.5+h/2btan⁡θK,DNU=Wcos⁡θ0.5-h/2btan⁡θK,DS=Wsin⁡θK,where DNL and DNU are the normal embedding distances of the lower contact point PL and upper contact point PU, respectively, DS is the sum of all the tangential embedding distances, W is the weight of the block, h and b are the height and width of the block, respectively, and θ is the dip angle of the slope.Comparing the analytical resultd and DDA result d′, r=d-d′/d is taken as the relative error; the results are shown in Table 3.Table 3
Comparison of the embedding distances for block stabilized using the three calculation solutions.
Contact point Analytical solution Conventional DDA method solution Improved DDA method solution Embedding distance (10−4 m) Embedding distance (10−4 m) Relative error (%) Embedding distance (10−4 m) Relative error (%) DNL 2.0944 2.0759 0.8833 2.0684 1.1214 DNU 0.5304 0.5254 0.9426 0.5236 1.2255 DSL — 1.8186 — 1.8097 — DSU — 1.3094 — 1.3086 — DS 3.1280 3.1016 0.8439 3.0894 0.9463As shown in Table3, the difference between the DDA results and analytical results is small and the relative error is controlled within 1.5%, which illustrates the proposed method’s high accuracy. Thus, it can be assumed that the calculating accuracy of the contact force is high and can satisfy the engineering calculation requirements.


## 7.2. Model of a Tunnel

### 7.2.1. Numerical Model

The model of an underground tunnel with fully developed joints is shown in Figure12; the rock mass is cut into blocks by the joint planes. The blocks have the following material properties: bulk density γ=27KN/m3, Young’s modulus E=1.50 MPa, Poisson’s ratio μ=0.23, friction angle φ=5°, cohesive force c=0, and contact stiffness P=50E. All blocks are subjected to the same initial stresses of σx=-0.05 MPa, σy=-1.05 MPa, and τxy=0, and there is no body force.Figure 12
Model of the tunnel before deformation.


### 7.2.2. Calculation Results

The time step is taken asΔ=0.01 s, and the deformation of the tunnel after 100 time steps is illustrated in Figure 13. The surrounding rock is in a state of self-stabilization, and the deformation and displacement of the block in the four corners are large compared with those in other locations. The outer contour of the tunnel is fixed, and no displacement occurs.Figure 13
Model of the tunnel after deformation.The conventional DDA method and improved DDA method differ regarding the calculation process. The difference mainly concerns the contact detection and solving of the extremum of the potential energy, not the pretreatment, postprocessing, and analysis of the potential energy inside the blocks. In the conventional DDA method, implicit open-close iteration and solving the large sparse linear equations are required at every time step. Implicit open-close iteration can be viewed as the process of determining the contact constraint, whereas solving the large sparse linear equations can be seen as the process of analyzing the extremum of the potential energy. More time is consumed conducting implicit open-close iteration compared with solving large sparse linear equations. In contrast, in the improved DDA method, there is no need to conduct implicit open-close iteration, and the contact constraints between blocks have already been expressed using the approximated step function and Lagrange interpolation. The expression of the overall potential energy function of the improved DDA method is more complex than that of the conventional DDA method, and the time is mainly consumed while determining the extremum of the potential energy.Because the overall potential energy function and the process of solving equations differ between the two DDA methods, the iteration step cannot be used to evaluate the efficiency of the two methods. In this paper, the time consumed during calculation is used as the comparison criterion. Because the pretreatment, postprocessing, and analysis of the potential energy inside blocks of the two DDA methods are the same and their difference with respect to time consumed is mainly caused by contact detection and solving the extremum of the potential energy, it is reasonable to take the time consumed as the criterion. All analyses are performed on the same computer with the following system configuration: processor speed = 3.40 GHz and RAM = 8.0 GB.Four calculation conditions, with contact stiffness of 80E, 100E, 120E, and 150E, are given in Table 4. The conventional DDA method and improved DDA method are used to analyze the stability of the surrounding rock of the tunnel after 100 time steps, and the time step is taken as Δ=0.01 s. The time consumed is recorded, and the results are shown in Table 4.Table 4
Comparison of iteration time with different contact stiffness.
Contact stiffness Iteration time (s) Conventional DDA method Improved DDA method 5 0 E 793.2 592.6 80 E 674.9 602.4 100 E 812.7 613.7 120 E Terminated 589.5 150 E 931.9 598.3As shown in Table4, the time consumed and the variation trends of the two DDA methods are different. The improved DDA method is more stable, and the time consumed for different contact stiffness remains at approximately 600 s; the contact stiffness has only a slight impact on the time consumed. Compared with the improved DDA method, the time consumed using the conventional DDA method is generally over 600 s, the variation trends are not distinct, and the embedding distance varies with changing contact stiffness. Because the contact constraint is uncertain in the conventional DDA method, considerable time is spent performing the implicit open-close iteration. If the open-close iteration is successful, it needs to be conducted only a few times to obtain the precise contact constraints; thus, the amount of time consumed is not overly burdensome. If the precise contact constraint cannot be obtained after conducting the implicit open-close iteration more than six times, the time step will be reduced to 1/3 of the original. In this case, the time step will be smaller and considerable time will be spent in the previous process; thus, considerably more time is consumed by the conventional DDA method than by the improved DDA method.When the value of the contact stiffness is 120E, the calculation process with the conventional DDA method cannot be completed. The sign “terminated” means that, in the process of implicit open-close iteration, the precise contact constraints could not be found and the time step was reduced several times, finally falling into an infinite loop and being forced to terminate [33]. The robustness of the conventional DDA method is not strong, and in some cases, the implicit open-close iteration cannot be conducted successfully, thus preventing the final correct solution from being obtained. Because implicit open-close iteration is not needed in the improved DDA method, the above problems do not occur. Thus, the improved DDA method has a more stable and more robust calculation convergence.Based on the two aspects above, the improved DDA method exhibits a faster calculation speed than the conventional DDA method; moreover, it has a more robust calculation convergence.


## 7.2.1. Numerical Model

The model of an underground tunnel with fully developed joints is shown in Figure12; the rock mass is cut into blocks by the joint planes. The blocks have the following material properties: bulk density γ=27KN/m3, Young’s modulus E=1.50 MPa, Poisson’s ratio μ=0.23, friction angle φ=5°, cohesive force c=0, and contact stiffness P=50E. All blocks are subjected to the same initial stresses of σx=-0.05 MPa, σy=-1.05 MPa, and τxy=0, and there is no body force.Figure 12
Model of the tunnel before deformation.


## 7.2.2. Calculation Results

The time step is taken asΔ=0.01 s, and the deformation of the tunnel after 100 time steps is illustrated in Figure 13. The surrounding rock is in a state of self-stabilization, and the deformation and displacement of the block in the four corners are large compared with those in other locations. The outer contour of the tunnel is fixed, and no displacement occurs.Figure 13
Model of the tunnel after deformation.The conventional DDA method and improved DDA method differ regarding the calculation process. The difference mainly concerns the contact detection and solving of the extremum of the potential energy, not the pretreatment, postprocessing, and analysis of the potential energy inside the blocks. In the conventional DDA method, implicit open-close iteration and solving the large sparse linear equations are required at every time step. Implicit open-close iteration can be viewed as the process of determining the contact constraint, whereas solving the large sparse linear equations can be seen as the process of analyzing the extremum of the potential energy. More time is consumed conducting implicit open-close iteration compared with solving large sparse linear equations. In contrast, in the improved DDA method, there is no need to conduct implicit open-close iteration, and the contact constraints between blocks have already been expressed using the approximated step function and Lagrange interpolation. The expression of the overall potential energy function of the improved DDA method is more complex than that of the conventional DDA method, and the time is mainly consumed while determining the extremum of the potential energy.Because the overall potential energy function and the process of solving equations differ between the two DDA methods, the iteration step cannot be used to evaluate the efficiency of the two methods. In this paper, the time consumed during calculation is used as the comparison criterion. Because the pretreatment, postprocessing, and analysis of the potential energy inside blocks of the two DDA methods are the same and their difference with respect to time consumed is mainly caused by contact detection and solving the extremum of the potential energy, it is reasonable to take the time consumed as the criterion. All analyses are performed on the same computer with the following system configuration: processor speed = 3.40 GHz and RAM = 8.0 GB.Four calculation conditions, with contact stiffness of 80E, 100E, 120E, and 150E, are given in Table 4. The conventional DDA method and improved DDA method are used to analyze the stability of the surrounding rock of the tunnel after 100 time steps, and the time step is taken as Δ=0.01 s. The time consumed is recorded, and the results are shown in Table 4.Table 4
Comparison of iteration time with different contact stiffness.
Contact stiffness Iteration time (s) Conventional DDA method Improved DDA method 5 0 E 793.2 592.6 80 E 674.9 602.4 100 E 812.7 613.7 120 E Terminated 589.5 150 E 931.9 598.3As shown in Table4, the time consumed and the variation trends of the two DDA methods are different. The improved DDA method is more stable, and the time consumed for different contact stiffness remains at approximately 600 s; the contact stiffness has only a slight impact on the time consumed. Compared with the improved DDA method, the time consumed using the conventional DDA method is generally over 600 s, the variation trends are not distinct, and the embedding distance varies with changing contact stiffness. Because the contact constraint is uncertain in the conventional DDA method, considerable time is spent performing the implicit open-close iteration. If the open-close iteration is successful, it needs to be conducted only a few times to obtain the precise contact constraints; thus, the amount of time consumed is not overly burdensome. If the precise contact constraint cannot be obtained after conducting the implicit open-close iteration more than six times, the time step will be reduced to 1/3 of the original. In this case, the time step will be smaller and considerable time will be spent in the previous process; thus, considerably more time is consumed by the conventional DDA method than by the improved DDA method.When the value of the contact stiffness is 120E, the calculation process with the conventional DDA method cannot be completed. The sign “terminated” means that, in the process of implicit open-close iteration, the precise contact constraints could not be found and the time step was reduced several times, finally falling into an infinite loop and being forced to terminate [33]. The robustness of the conventional DDA method is not strong, and in some cases, the implicit open-close iteration cannot be conducted successfully, thus preventing the final correct solution from being obtained. Because implicit open-close iteration is not needed in the improved DDA method, the above problems do not occur. Thus, the improved DDA method has a more stable and more robust calculation convergence.Based on the two aspects above, the improved DDA method exhibits a faster calculation speed than the conventional DDA method; moreover, it has a more robust calculation convergence.


## 8. Conclusions

This paper proposed an improved DDA method based on explicitly solving contact constraints to solve the problem of poor convergence in implicit open-close iteration; several conclusions can be drawn from this study:(1)
The approximated step functions are used to represent the relationship between contact constraints and displacement variables. The overall potential energy function containing only displacement variables as an unknown is deduced based on Lagrange interpolation. After determining the block displacement using the variable metric method, the current contact constraints can be explicitly obtained by substituting the displacement variables into the approximated step function, which can help avoid using implicit open-close iteration.(2)
The numerical examples illustrate that the error range and variation trends of the conventional DDA method and improved DDA method are highly similar. The improved DDA method can control the calculation accuracy with the approximation factorε and termination condition ∇ΠDk≤ε′. Thus, it can be assumed that the improved DDA method has a high calculation accuracy. In contrast, the conventional DDA method requires implicit open-close iteration, which has poor convergence; because there is no need to consider implicit open-close iteration in the improved DDA method, it has a more robust calculation convergence. Thus, the improved DDA method based on explicitly solving contact constraints is a more stable and effective numerical method for solving the problem of discontinuity and provides a new approach to improving open-close iteration.



---

*Source: 9452175-2016-10-17.xml*