Improved DDA Method Based on Explicitly Solving Contact Constraints

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.


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 2 Mathematical Problems in Engineering 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 openclose 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][17][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-toedge 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.

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 point   1.
Embedding can be expressed as Equation (1) illustrates that Δ is a vector-valued function containing only the displacement variables {} as an unknown, which can be expressed as Δ = Δ ({}) . ( If Δ is smaller than 0, embedding occurs and normal and tangential embedding distances exist.The normal embedding distance   and the tangential embedding distance   can be expressed as [1] where {  } and {  } are the displacement variables of the two blocks in contact with each other and { 1 }, { 2 }, { 1 }, { 2 },  1 ,  2 , and  are the intermediate variables. Equation (3) illustrates that the normal embedding distance   and the tangential embedding distance   are vectorvalued function containing only the displacement {} as an unknown, which can be expressed as Variable Δ, normal embedding distance   , and tangential embedding distance   are deduced for the following deduction of the approximated step function.

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 Table 1.
As shown in Table 1, 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 (  < 0) and 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 (  < 0) and ( ⃗  ⋅ ⃗  < 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.

Approximated
Step Function  1 (Δ).As shown in Table 1, when   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  1 (Δ) 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  1 (Δ) can be expressed as where Δ reflects the degree of embedding and is expressed as (1).The step function  1 (Δ) is illustrated in Figure 2. The specific value of Δ is not known before the calculation; thus, the specific value of  1 (Δ) 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  1 (Δ); they are able to fully approximate the step function, whose value ranges from −1 to 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.
In Figure 3, the horizontal coordinate Δ reflects the degree of embedding, and the vertical coordinate (Δ) 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.
Figure 3 illustrates that the degree of approximation of hyperbolic tangent function  2 (Δ) = tanh (Δ/) is high.When Δ ranges from 0 to 10 −7 , the slope of the curve of function  2 is high and the value grows rapidly.When Δ reaches 10 −7 , the curve of function  2 remains fairly flat.From this aspect, the degree of approximation of function  2 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,  1 (Δ),  3 (Δ), and  4 (Δ), are rather low compared to that of  2 (Δ).The expression of the hyperbolic tangent function  2 (Δ) is relatively simple to express and easy to calculate.Thus, the hyperbolic tangent function  2 (Δ) is used to approximate the step function  1 (Δ), which can be written as To ensure that (6) can fully approximate the step function  1 (Δ), the approximation factor  in ( 6) should be sufficiently small and the magnitude should approach 10 −8 .Thus,  1 (Δ) can be used to simulate whether the current type of contact constraint is open or embedding.

Approximated
Step Function  2 (  /  ).Similarly, embedding is known to occur in the current time step when   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 function  2 (  /  ) can be used to simulate the fixed embedding and sliding embedding, in which   is the tangential embedding distance and   is the normal embedding distance.When calculating tangential embedding distance, we cannot determine the sign of   , only the absolute value.Suppose that when   /  is larger than tan  or smaller than − tan , the value of the function is 1 and sliding embedding occurs; when   /  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  2 (  /  ) is illustrated in Figure 4 and it can be expressed as With reference to (6),  2 (  /  ) can be approximated using the following equation: Thus,  2 (  /  ) 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.

Approximated
Step Function  3 ( ⃗  ⋅ ⃗ ).When the type of contact of the last time step is sliding, another step function Sliding embedding  3 ( ⃗  ⋅ ⃗ ) can be used to simulate the fixed embedding and sliding embedding, in which ⃗  is the friction vector and ⃗  is the shear displacement vector.The friction vector ⃗  is the product of the tangential stiffness  and the shear displacement vector  →   of the last time step, and the shear displacement vector ⃗  of the current time step is the vector of the tangential embedding distance   in Section 2. Suppose that when ⃗  ⋅ ⃗  is larger than 0, the value of the function is 1 and sliding embedding occurs, whereas when   /  is smaller than 0, the value of the function is −1 and fixed embedding occurs.The function  3 ( ⃗  ⋅ ⃗ ) is illustrated in Figure 5 and it can be expressed as With reference to (6),  3 ( ⃗  ⋅ ⃗ ) can be approximated using the following equation: Thus,  3 ( ⃗  ⋅ ⃗ ) 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 Δ,   ,   , and ⃗  are functions containing only displacement variables {} as an unknown; thus, ( 6), (8), and (10) are also functions containing only displacement variables {} 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 {} = ({}) can be generalized using  1 (Δ),  2 (  /  ), and  3 ( ⃗ ⋅ ⃗ ).Next, the concept of Lagrange interpolation is introduced to deduce the expression of the contact potential energy function Π c .

Contact Potential Energy Function Based on Lagrange Interpolation
Approximated step functions  1 (Δ),  2 (  /  ), and  3 ( ⃗ ⋅ ⃗ ) 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.
The potential energy function in response to different types of contact constraints is not the same but has a oneto-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 When the type of contact of the last time step is sliding, it can be written as Mathematical Problems in Engineering Equations ( 11) and ( 12) represent the current contact potential energy function with different types of contact constraints of the last time step, where  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 {} as an unknown.
Suppose that the number of contact pairs is .Then, (13) can be deduced by summing up all the contact potential energy values.
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 The contact constraints have already been represented using the displacement of blocks, and ( 14) contains only the displacement variables {} as an unknown 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.

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.

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: where { (+1) } and { () } are the displacement variables of blocks at time steps  + 1 and , respectively, Π is the overall potential energy, ∇Π is the gradient of the overall potential energy with respect to the displacement variables of blocks {},  () is the optimal step length of time step , and  () is a symmetric matrix with size  ×  (where  is the number of freedom of the overall displacement variables).Next, the iteration factors mentioned above are discussed in detail.

Determination of Iteration Factors.
The three main iteration factors are the gradient of the overall potential energy ∇Π, the optimal step length factor  () , and the symmetric matrix  () .Currently, the ahead difference quotient is widely used to approximate the gradient, with ∇Π({ () }) replacing ∇Π({ () }); the expression can be written as where ℎ  controls the accuracy of the calculation and  is the unit vector in -dimensional Euclidean space. () is the optimal step length factor of time step .The basic process of realizing this method is determining the value of the optimal step length  () 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 matrix  () can fully approximate the inverse matrix of the Hesse matrix [∇ 2 Π({ () })] −1 when reaching the extremum.The variable metric matrix of the adjacent time step has the following relationship: where Δ () is the calibrated matrix and can be obtained using the DFP method, written as follows: where { () } is the displacement variables of blocks at time step , {Δ () } is the displacement difference {Δ () } = { (+1) } − { () },  () is the gradient of the overall potential energy of the displacement of blocks at time step , and Δ ()  is the gradient difference Δ () =  (+1) −  () .

Process of Iteration.
The optimal solution of displacement variables is obtained with the following steps: (1) Determine the overall potential energy function ∇Π({}) deduced in Sections 3 and 4.

Open-close iteration
Open-close iteration not needed Time step i j = 1 j = 2 j = n, achieve the state of convergence Time step i + 1

Implementation Procedure
The difference in contact detection between the conventional DDA method and improved DDA method is illustrated in Figure 6.
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 Figure 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.

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 Figure 8; 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  1 is 8 m.The dip angle  2 of the lower part of the slope is 20 ∘ , and the length  2 is 18 m.Young's modulus  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.
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.

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   and the two types of DDA results  are analyzed. = (−  )/ is taken as the relative error.The results are shown in Figures 9-11.
As shown in Figure 9, 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 Figure 9 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. Conventional DDA method Improved DDA method 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 6.0 6.5 0.0 Relative error between two DDA solutions (%) The relative error in Figure 10 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.0s.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 Figures 9, 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 ‖∇Π( () )‖ ≤   .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: where DN L and DN U are the normal embedding distances of the lower contact point  L and upper contact point  U , respectively, DS is the sum of all the tangential embedding distances,  is the weight of the block, ℎ and  are the height and width of the block, respectively, and  is the dip angle of the slope.
Comparing the analytical result  and DDA result   ,  = ( −   )/ is taken as the relative error; the results are shown in Table 3.
As shown in Table 3, 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.

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.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 openclose 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  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.
As shown in Table 4, 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 openclose 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.

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 openclose iteration. ( 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 ‖∇Π( () )‖ ≤   .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.

Figure 1 :
Figure 1: Vertex-edge contact in the block.

1 Figure 2 :
Figure 2: Schematic diagram of embedding simulated by a step function.

Figure 3 :
Figure 3: Approximated function curve of the step function.

Figure 4 :
Figure 4: Schematic diagram of the embedding mode simulated by a step function (open or fixed in the last calculation step).

Figure 5 :
Figure 5: Schematic diagram of the embedding mode simulated by a step function (sliding in the last calculation step).

Figure 6 :
Figure 6: Comparison of contact detection between the two DDA methods.

Figure 7 :Figure 8 :
Figure 7: Calculation process of the improved DDA method.

Figure 9 :
Figure 9: Absolute error-time curve between two DDA solutions and analytical solution.

Figure 10 :
Figure 10: Relative error-time curve between two DDA solutions and analytical solution (overview).

Figure 11 :
Figure 11: Relative error-time curve between two DDA solutions and analytical solution (partially zoomed-in view).

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  = 27 KN/m 3 , Young's modulus  = 1.50 MPa, Poisson's ratio  = 0.23, friction angle  = 5 ∘ , cohesive force  = 0, and contact stiffness  = 50.All blocks are subjected to the same initial stresses of   = −0.05MPa,   = −1.05MPa, and   = 0, and there is no body force.

Figure 13 :
Figure 13: Model of the tunnel after deformation.

Table 1 :
Discriminant table of contact constraints.

Table 2 :
Discriminant table of contact constraints with an approximated step function.

Table 4 :
Comparison of iteration time with different contact stiffness.