Improved DDA Method Based on Explicitly Solving Contact Constraints

Zhao, Jian; Xiao, Ming; Chen, Juntao; Li, Dongdong

doi:https://doi.org/10.1155/2016/9452175

Mathematical Problems in Engineering

On this page

Abstract Introduction Conclusions Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2016 | Article ID 9452175 | https://doi.org/10.1155/2016/9452175

Improved DDA Method Based on Explicitly Solving Contact Constraints

Jian Zhao,^1,2Ming Xiao ,^1,2Juntao Chen,^1,2and Dongdong Li^1,2

Academic Editor: Alberto Borboni

Received11 Jul 2016

Accepted25 Aug 2016

Published17 Oct 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.

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 point moves to the position of point after deforming, is the contact line, and and are the coordinate values and displacement of point , respectively. Embedding occurs if points , , and rotate from to in the same direction and point crosses over the line , which is illustrated in Figure 1.

(a) Before embedding

(b) After embedding

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 , , , , , , and are the intermediate variables.

Equation (3) illustrates that the normal embedding distance and the tangential embedding distance are vector-valued 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.

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 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 and 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 and 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 Function

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 simulating whether embedding occurs: when is smaller than 0, the value of the function is and embedding does not occur; when is larger than 0, the value of the function is and the embedding occurs. The step function can be expressed aswhere reflects the degree of embedding and is expressed as (1). The step function is illustrated in Figure 2.

The specific value of is not known before the calculation; thus, the specific value of 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 ; they are able to fully approximate the step function, whose value ranges from to :(1)Square root function (2)Hyperbolic tangent function (3)Arctangent function (4)Exponential function

is the approximation factor. When approaches , the degree of approximation is extremely high. The curves of the three approximated step functions, with being , 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⁻⁷, 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 is high. When ranges from 0 to 10⁻⁷, the slope of the curve of function is high and the value grows rapidly. When reaches 10⁻⁷, the curve of function remains fairly flat. From this aspect, the degree of approximation of function 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, , , and , are rather low compared to that of . The expression of the hyperbolic tangent function is relatively simple to express and easy to calculate. Thus, the hyperbolic tangent function is used to approximate the step function , which can be written asTo ensure that (6) can fully approximate the step function , the approximation factor in (6) should be sufficiently small and the magnitude should approach 10⁻⁸.

Thus, can be used to simulate whether the current type of contact constraint is open or embedding.

3.2. Approximated Step Function

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 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 or smaller than , the value of the function is and sliding embedding occurs; when is larger than and smaller than , the value of the function is and fixed embedding occurs, where is the friction angle. The function is illustrated in Figure 4 and it can be expressed asWith reference to (6), can be approximated using the following equation:

Thus, 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 Function

When the type of contact of the last time step is sliding, another step function 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 , the value of the function is and sliding embedding occurs, whereas when is smaller than , the value of the function is and fixed embedding occurs. The function is illustrated in Figure 5 and it can be expressed asWith reference to (6), can be approximated using the following equation:

Thus, 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 , , and . Next, the concept of Lagrange interpolation is introduced to deduce the expression of the contact potential energy function .

4. Contact Potential Energy Function Based on Lagrange Interpolation

Approximated step functions , , and can yield a certain value ( to ), 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 one-to-one correspondence. When fixed embedding occurs, there exists a fixed contact potential energy function ; when sliding embedding occurs, there exists a sliding contact potential energy function . 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 asWhen the type of contact of the last time step is sliding, it can be written as

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 and sliding contact potential energy 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 of the block system. The potential energy inside blocks is generated with stress, strain, and loading. For details regarding the block potential energy , refer to the conventional DDA method [1]. The overall potential energy can be deduced by summing up the contact potential energy and the block potential energy 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.

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:where and are the displacement variables of blocks at time steps 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.

5.2. 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 aswhere 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 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 , is the gradient of the overall potential energy of the displacement of blocks at time step , and is the gradient difference .

5.3. 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.(2)Specify the initial iteration point and termination factor .(3)Set time step being 0 and being identity matrix .(4)Calculate the value of referred to in (17) and set being .(5)Set the value of in order to minimize the value of which can be referred to in [29].(6)Check if the termination condition is verified. If so, stop iterating and set being the optimal solution of displacement variables. If not, conduct the following step.(7)Calculate the value of and referred to in (18) and (19). Set time step being and turn to step .

After we obtain the displacement of blocks , the contact constraints can be obtained by substituting into . The relationship between displacement variables and contact constraints between blocks has been generalized using , , and 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 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.

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 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 of the upper part of the slope is 50°, and the length is 8 m. The dip angle of the lower part of the slope is 20°, and the length 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.

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 but larger than the dip angle of the lower part ; 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.

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.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 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 and are the normal embedding distances of the lower contact point and upper contact point , respectively, 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.

7.2. Model of a Tunnel

7.2.1. Numerical Model

The model of an underground tunnel with fully developed joints is shown in Figure 12; the rock mass is cut into blocks by the joint planes. The blocks have the following material properties: bulk density , Young’s modulus MPa, Poisson’s ratio , friction angle , cohesive force , and contact stiffness . All blocks are subjected to the same initial stresses of MPa, MPa, and , and there is no body force.

7.2.2. Calculation Results

The time step is taken as 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 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 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 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 . 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.

Competing Interests

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

Acknowledgments

This work was supported by the National Key Basic Research and Development Plan (973) of China (Project no. 2015CB057904) and the National Natural Science Foundation of China (Project no. 51579191). This support is acknowledged and greatly appreciated.

References

G. H. Shi, Discontinuous deformation analysis—a new numerical model for the statics and dynamics of block systems [Ph.D. thesis], Department of Civil Engineering, University of California Berkeley, Berkeley, Calif, USA, 1988.
C. Y. Koo and J. C. Chern, “The development of DDA with third order displacement function,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, pp. 342–350, Albuquerque, NM, USA, 1996.
View at: Google Scholar
Y. M. Max, “Discontinuous deformation analysis the third order displacement function,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA '96) and Simulations of Discontinuous Media, pp. 383–395, Albuquerque, NM, USA, 1996.
View at: Google Scholar
K. Shyu and M. R. Salami, “Stability of bartlett dam,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, pp. 440–446, Albuquerque, NM, USA, 1996.
View at: Google Scholar
W. C. Yue, C. C. Chuang, and T. S. Jopan, “Time intergration theories for the DDA method with finite element meshes,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, pp. 263–288, Albuquerque, NM, USA, 1996.
View at: Google Scholar
B. Amadei, L. Chihsen, and J. Dwyer, “Recent extensions to the DDA method,” in Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, pp. 1–30, Albuquerque, NM, USA, 1996.
View at: Google Scholar
S. A. R. Beyabanaki, B. Ferdosi, and S. Mohammadi, “Validation of dynamic block displacement analysis and modification of edge-to-edge contact constraints in 3-D DDA,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 7, pp. 1223–1234, 2009.
View at: Publisher Site | Google Scholar
C. T. Lin, B. Amadei, J. Jung, and J. Dwyer, “Extensions of discontinuous deformation analysis for jointed rock masses,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 33, no. 7, pp. 671–694, 1996.
View at: Publisher Site | Google Scholar
C. Y. Koo and J. C. Chern, “Modification of the DDA method for rigid block Problems,” International Journal of Rock Mechanics and Mining Sciences, vol. 35, no. 6, pp. 683–693, 1998.
View at: Publisher Site | Google Scholar
M. M. MacLaughlin and D. M. Doolin, “Review of validation of the discontinuous deformation analysis (DDA) method,” International Journal for Numerical & Analytical Methods in Geomechanics, vol. 30, no. 4, pp. 271–305, 2006.
View at: Publisher Site | Google Scholar
A. Borboni, M. Mor, and R. Faglia, “Gloreha—hand robotic rehabilitation: design, mechanical model, and experiments,” Journal of Dynamic Systems, Measurement, and Control, vol. 138, no. 11, Article ID 111003, 12 pages, 2016.
View at: Publisher Site | Google Scholar
Q. Tian, Z. Zhao, and H. Bao, “Block fracturing analysis using nodal-based discontinuous deformation analysis with the double minimization procedure,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 9, pp. 881–902, 2014.
View at: Publisher Site | Google Scholar
Y. M. Cheng, “Advances and improvements in discontinuous deformation analysis,” Computers and Geotechnics, vol. 22, no. 1, pp. l53–I63, 1998.
View at: Publisher Site | Google Scholar
B. Amadei, C. Lin, S. Sture et al., “Modeling fracture of rock masses with the DDA method,” in Rock Mechanics, L. Nelson, Ed., pp. 583–590, Balkema, Rotterdam, The Netherlands, 1994.
View at: Google Scholar
H. Bao and Z. Zhao, “The vertex-to-vertex contact analysis in the two-dimensional discontinuous deformation analysis,” Advances in Engineering Software, vol. 45, no. 1, pp. 1–10, 2012.
View at: Publisher Site | Google Scholar
S. A. R. Beyabanaki, R. G. Mikola, and K. Hatami, “Three-dimensional discontinuous deformation analysis (3-D DDA) using a new contact resolution algorithm,” Computers and Geotechnics, vol. 35, no. 3, pp. 346–356, 2008.
View at: Publisher Site | Google Scholar
L. He, X. M. An, and Z. Y. Zhao, “Development of contact algorithm for three-dimensional numerical manifold method,” International Journal for Numerical Methods in Engineering, vol. 97, no. 6, pp. 423–453, 2014.
View at: Publisher Site | Google Scholar
J. Wu, “New edge-to-edge contact calculating algorithm in three-dimensional discrete numerical analysis,” Advances in Engineering Software, vol. 39, no. 1, pp. 15–24, 2008.
View at: Publisher Site | Google Scholar
H. Zheng and X. Li, “Mixed linear complementarity formulation of discontinuous deformation analysis,” International Journal of Rock Mechanics and Mining Sciences, vol. 75, pp. 23–32, 2015.
View at: Publisher Site | Google Scholar
H. Zheng and W. Jiang, “Discontinuous deformation analysis based on complementary theory,” Science in China, Series E: Technological Sciences, vol. 52, no. 9, pp. 2547–2554, 2009.
View at: Publisher Site | Google Scholar
X. X. Feng, W. U. Ai-Qing, and L. U. Bo, “Open-close iterations for DDA contacts,” Journal of Yangtze River Scientific Research Institute, vol. 31, no. 4, pp. 52–58, 2014.
View at: Google Scholar
A.-Q. Wu, X.-X. Feng, and B. Lu, “Parametric research on time step and spring stiffness in DDA,” Rock and Soil Mechanics, vol. 36, no. 3, pp. 891–897, 2015.
View at: Publisher Site | Google Scholar
Y. Zhang and Y. Cheng, “The extension and application of DDA method,” Chinese Journal of Geotechnical Engineering, vol. 20, no. 2, pp. 109–111, 1998.
View at: Google Scholar
P. A. Cundall, “Formulation of a three-dimensional distinct element model—part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 25, no. 3, pp. 107–116, 1988.
View at: Publisher Site | Google Scholar
B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” Journal of the Optical Society of America A: Optics and Image Science, and Vision, vol. 11, no. 4, pp. 1491–1499, 1994.
View at: Publisher Site | Google Scholar
K. Bathe and P. Bouzinov, “On the constraint function method for contact problems,” Computers & Structures, vol. 64, no. 5-6, pp. 1069–1085, 1997.
View at: Publisher Site | Google Scholar
P. Wriggers and T. A. Laursen, Computational Contact Mechanics, John Wiley & Sons, New York, NY, USA, 2003.
J. T. Chen, M. Xiao, and G. Liu, “Application of improved variable-metric method in displacement back analysis of underground engineering,” Chinese Journal of Underground Space and Engineering, vol. 4, no. 5, pp. 825–829, 2008.
View at: Google Scholar
E. K. P. Chong and S. H. Zak, Single-Dimensional Search Methods. Wiley Encyclopedia of Operations Research and Management Science, John Wiley & Sons, New York, NY, USA, 2011.
A. Wu, X. Ding, B. Lu, Q. Zhang, and G. Shi, “Validation for rock block stability and its application to rock slope stability evaluation using DDA method,” Chinese Journal of Rock Mechanics and Engineering, vol. 27, no. 4, pp. 664–672, 2008.
View at: Google Scholar
L. I. Yong, X. Feng, M. Luan, and Y. J. Wang, “Contact force element of discontinuous deformation computational mechanics model,” Journal of Northeastern University (Natural Science), vol. 23, pp. 1112–1115, 2002.
View at: Google Scholar
M. R. Yeung, Application of Shi's discontinuous deformation analysis to the study of rock behavior [Ph.D. thesis], University of California, Berkeley, Calif, USA, 1991.
W. Jiang, Discontinuous Deformation Analysis Method Based on Complementary Theory, Institute of Rock and Soil Mechanics, Graduate School of Chinese Academy of Sciences, 2009.

Copyright

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

PDF Download Citation

Download other formats

Order printed copies

Views

1124

Downloads

606

Citations