Explicit Dynamic DDA Method considering Dynamic Contact Force

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

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

Shock and Vibration

On this page

Abstract Introduction Discussion Conclusion Acknowledgments References Copyright Related Articles

Research Article | Open Access

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

Explicit Dynamic DDA Method considering Dynamic Contact Force

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

Academic Editor: Vadim V. Silberschmidt

Received16 Sept 2016

Revised22 Nov 2016

Accepted28 Nov 2016

Published28 Dec 2016

Abstract

This paper proposes an explicit dynamic DDA method considering dynamic contact force, which aims at solving the problems of low efficiency of dynamic contact detection and the simulation of dynamic contact force in the conventional DDA method. The mutual contact between blocks can be regarded as the application of point loading on a single block, and the corresponding contact submatrix can be calculated and the simultaneous equations of the block system can be integrated. The central difference method is adopted to deduce the explicit expression of block displacement containing dynamic contact force. With the relationship between displacement and dynamic contact force, contact constraint equations of a block system are obtained to calculate the dynamic contact force and the corresponding block displacement. The accuracy of the explicit dynamic DDA method is verified using two numerical cases. The calculation results show that the new DDA method can be applied in large-scale geotechnical engineering.

1. Introduction

The discontinuous deformation analysis method proposed by Shi [1] is a new and efficient numerical method for discontinuous rock mass. The DDA method can efficiently simulate mutual contact between blocks, accounting for the impact of geological joints on block movement. The DDA method is widely applied in geotechnical engineering.

Currently, dynamic issues such as impact, explosion, and earthquake are very common in geotechnical engineering. Compared with static loading, dynamic loading often changes over time. The inertial force of structure caused by acceleration cannot be ignored compared with the loading itself, and the degree of deformation caused by dynamic loading is often larger and the failure mode is often more complex. Much research work focuses on the study of dynamic problems using the DDA method. Pei et al. [2] simulated the process of a bridge pier’s dynamic response after being hit by an instable rock mass; Ma et al. [3] simulated the failure process of a Brazilian disc under dynamic loading; Ning et al. [4] simulated the phenomenon of a blast hole’s expansion, the failure of a rock mass, and the formation of a rock pile; Zhang et al. [5] simulated the law of propagation and attenuation of stress wave caused by blasting in a jointed rock mass; Ning [6, 7] simulated the failure process of a continuous medium with artificial joints and applied it to blasting engineering; Wu [8, 9] simulated the failure process of the Chiu-fen-erh-shan slope and the Tsaoling slope under seismic loading; Shi [10] studied the rockfall and collapse of the Yubabuna tunnel under seismic loading; Kong et al. [11] studied the stability of a concrete-faced rockfill dam under seismic loading. Moreover, the DDA method has been applied in the dynamic analysis of a bridge [12], retaining wall [13], towel [14], and levee [15].

Dynamic contact [16, 17] is one of the main issues studied in the dynamic DDA method, and the research focuses on the dynamic contact detection and the calculation of dynamic contact force. From the aspect of contact detection, the contact modes between blocks change rapidly under dynamic loading and contact and separation can occur repeatedly. A large number of contact predictions and judgments are needed and an open-close iteration is adopted in the conventional DDA method to distinguish the contact modes (fixed, sliding, and separation modes). In the process of open-close iteration, the contact modes are assumed in advance and corrected for those contact pairs that violate the contact criteria. Obviously, open-close iteration does not have strict mathematical convergence criteria, and the scheme cannot guarantee that the iteration is always convergent [18]. Meanwhile, the time consumed during the open-close iteration is very long. If the number of time steps of the dynamic analysis is very large, open-close iteration would consume much time and decrease the efficiency of calculation.

From the aspect of dynamic contact force, short embedding distance is allowed in the conventional DDA method, and contact force is the product of embedding distance and contact stiffness. The determination of the contact stiffness value heavily relies on the engineer’s experience and does not have a specific rule [19]. Meanwhile, under dynamic loading, the dynamic contact force between blocks is related to the impact velocity, the stress wave velocity, and Young’s modulus [20]. The product of embedding distance and contact stiffness in the conventional DDA method cannot reflect the dynamic feature. If the dynamic contact force cannot be calculated precisely, the corresponding block displacement cannot be accurate either.

The current numerical calculation methods aiming at the dynamic contact problems are the penalty method [21, 22] and the Lagrange multiplier method [23, 24] and their improved versions [25, 26]. The product of embedding distance and contact stiffness is regarded as the contact force in the penalty method, and it cannot represent the precise value of the dynamic contact force. The contact force is set as the unknown in the Lagrange multiplier method to meet the contact conditions precisely, but it increases the degree of freedom of the simultaneous equations and the difficulty of solving the simultaneous equations. Meanwhile, all of these methods adopt the open-close iteration to obtain the contact modes of the contact pairs, which consumes considerable time and decreases the efficiency of calculation. Aiming at the dynamic response of cracks, Liu [27, 28] proposed the dynamic contact force method, which has been applied in geotechnical engineering. The explicit scheme is adopted in the algorithm, and the displacement at the next time step can be expressed with the displacement and the contact force at the current time step. The dynamic solution of displacement can be obtained with the relationship between displacement and contact force. Due to the explicit scheme, the algorithm can greatly reduce the calculation time and computer memory consumed. The algorithm can easily reach convergence and is suitable for analyzing large-scale, complex dynamic contact problems. However, the contact pairs of the algorithm are independent and have no connection with each other, which violates the contact conditions in the DDA method. The contact pairs in the block system in the DDA method can influence each other, and the modification of one contact pair’s mode may lead to the modification of another contact pair’s mode. From this perspective, the dynamic contact force method cannot precisely simulate the contact mode of the complex block system. Moreover, there are other numerical methods, such as the impulse model method [29], the initial displacement method [30], and the dynamic contact model method [31]. These methods are mainly applied in the analysis of a crack’s dynamic response and are rarely applied in geotechnical engineering.

This paper proposes the explicit dynamic DDA method considering dynamic contact force based on the previous research, which aims at solving the problems of low efficiency of dynamic contact detection and the simulation of dynamic contact force in the conventional DDA method. Because there are many factors that can affect the dynamic contact force and it is difficult to determine the calculation formula of the dynamic contact force, with the reference of Liu’s dynamic contact force method, dynamic contact force is set as an unknown in this paper. The explicit expression of block displacement containing dynamic contact force is deduced and the contact constraint equations of a block system are obtained based on the contact relationship between displacement and dynamic contact force. The simultaneous equations are solved to obtain the corresponding displacement and dynamic contact force. There is no need to conduct the open-close iteration repeatedly to obtain the contact modes, and the dynamic contact force can be precisely calculated. Due to the explicit scheme, it can reduce the calculation time and the computer memory consumed. The calculating efficiency is greatly improved, and it can be applied in large-scale geotechnical engineering containing complex dynamic contact problems.

2. The Simultaneous Equations of the Explicit Dynamic DDA Method considering Dynamic Contact Force

The conventional DDA method is based on the principle of minimum potential energy. The potential energy of the block system is calculated and integrated, and the displacement can be obtained by minimizing the overall potential energy, which is the same as for the finite element method. However, the inertia force potential energy generated by unbalanced force and the contact potential energy generated by the contact between blocks must be calculated in the conventional DDA method, while this is not needed in the finite element method. In this section, the simultaneous equations of the explicit dynamic DDA method are deduced from the two aspects of inertia force potential energy and contact potential energy.

2.1. Dynamic Contact Force Simulated by Point Loading

A block system is composed of many single blocks, and contact force exists between blocks. The overall simultaneous equations of a block system can be obtained by integrating the simultaneous equations of every single block. In the following deduction, contact force is set as an unknown. Assume that the block contacts the adjacent blocks; the contact effect between the block and the adjacent blocks can be regarded as the point loading, which is illustrated in Figure 1.

The contact force applied on block as point loading can be classified into two categories. The first category is active contact force , as illustrated in Figure 1, and it is applied to the contact point of block itself. The second category is passive contact force , as illustrated in Figure 1, and it is applied to the contact line of block itself. From the principle of interacting force, the passive contact force applied to block and the active contact force applied to block are a pair of interacting forces; thus, the relationship between the two interacting forces can be expressed as

The active contact force can be decomposed into in the direction of and in the direction of , and the contact potential energy generated by active contact force applied to block can be expressed as where is the displacement of block and is the displacement function at the contact point of block . Taking the derivative of on displacement , a submatrix of point loading can be obtained and expressed as

Meanwhile, the point loading submatrix corresponding to passive contact force can be expressed as To unify the context, all of the equations are deduced with active contact force. With reference to (1), (4) can be expressed as

After the contact submatrix of block is calculated, the remaining submatrices, including elastic submatrix, initial stress submatrix, loading submatrix, and displacement-constraint submatrix, except the inertia force submatrix of block , are calculated and integrated into the stiffness matrix and the loading matrix . With reference to (4) and (5), the simultaneous equations of block can be expressed as

Assume that the number of blocks of the whole block system is . A similar deduction can be made for every block, and the simultaneous equations can be expressed asIntegrating (7) of all blocks, the overall simultaneous equations of the block system can be expressed asAll of the contact force in (8) is taken as active contact force, and inertia force is not taken into consideration.

2.2. Explicit Expression of Block Displacement Containing Dynamic Contact Force

Substituting the inertia force submatrix into (8), the overall simultaneous equations of the block system considering inertia force can be expressed as where is the overall mass matrix of the block system after being integrated and the mass submatrix of each single block can be expressed as where is the density of a block and is the displacement function.

The acceleration at every time step can be assumed as a constant in the conventional DDA method and can be expressed as

The implicit scheme is adopted to solve (9) in the conventional DDA method, and (9) needs to be iteratively solved at every time step. Due to the large degree of freedom of the overall simultaneous equations, a large amount of disk space needs to be occupied and much time needs to be consumed during the process of calculation. The contact modes are unknown at the beginning of every time step, and open-close iteration needs to be conducted repeatedly to determine the contact modes. Thus, the overall simultaneous equations need to be solved repeatedly at every time step. Obviously, the implicit scheme leads to low calculation efficiency while solving the large-scale dynamic DDA simultaneous equations containing complex contact problems, which needs to be improved.

In this paper, instead of taking (11) to approximate the acceleration at every time step, the central difference method [32] is adopted. The central difference method is the explicit algorithm, and solving the equilibrium equations iteratively is not needed with the explicit algorithm. Its calculation speed is very fast and it demands less computation memory, which leads to high computational efficiency. The central difference method is very suitable for solving large-scale geotechnical problems with dynamic contact problems, and the acceleration can be expressed asSubstituting (12) into (9), the discretization of the overall simultaneous equations in the time domain can be expressed as The displacement at time step can be obtained using (13) and can be expressed as

It can be seen from (14) that the displacement at time step is determined by the displacement at time steps and , which is the distinctive feature of the explicit algorithm [32]. From (14), the displacement at time step consists of two parts [18]: the predictable displacement without considering dynamic contact force and the amended displacement caused by dynamic contact force. The predictable displacement is only related to the displacement and the external load at time steps and and is a known quantity. Thus, the dynamic contact force needs to be solved to obtain the displacement at time step .

3. Contact Constraint Equations of the Block System

In the conventional DDA method, open-close iteration is adopted to detect contact modes, and it is a test-and-error scheme. The scheme cannot guarantee that the iteration is always convergent. In this paper, after the explicit expression of block displacement containing dynamic contact force is deduced, the contact constraint between the dynamic contact force and the block displacement [33, 34] is introduced to form the contact constraint equations, which can play the role of open-close iteration.

3.1. Normal Contact Constraint

From the aspect of normal contact constraint, the contact mode of the contact pair may be open or embedding. The relationship between displacement and dynamic contact force can be illustrated as in Figure 2.

When the contact mode is embedding, (15) is met and expressed asWhen the contact mode is open, (16) is met and expressed as where is the normal relative distance and is the normal contact force. According to (15) and (16), (17) can be expressed as follows:

Normal relative distance can be expressed as where is the initial spacing of the contact pair at the beginning of the time step, and are the displacements of the two contacted blocks, respectively, and are the displacement matrices of the two nearest points of the contact pair at the beginning of the time step, and is the normal vector, which points towards the inner part of the block.

The contact force in (14) is presented in the global coordinate system, but the normal contact force in (17) is presented in the local coordinate system. Thus, the coordinate transform is needed between the two groups of contact force in (14) and (17), and the coordinate transform can be expressed as

3.2. Tangential Contact Constraint

From the aspect of tangential contact constraint, the contact mode of the contact pair may be sliding or fixed. The relationship between displacement and dynamic contact force can be illustrated as in Figure 3.

When the contact mode is sliding, (20) is met and expressed as When the contact mode is fixed, (21) is met and expressed as where is the tangential relative distance, is the normal contact force, is the tangential contact force, and is the friction coefficient. According to (20) and (21), (22) can be expressed as

Tangential relative distance can be expressed as where and are the displacements of the two contacted blocks, respectively, and are the displacement matrices of the two nearest points of the contact pair at the beginning of the time step, and is the tangential vector, which rotates from the normal vector clockwise.

Similarly, the coordinate transform is needed between the tangential contact force in the local coordinate system and the contact force in the global coordinate system. The coordinate transform can be expressed as

3.3. Smoothing Contact Constraint Equations

Each contact pair is jointly controlled by the normal contact constraint and tangential contact constraint. According to (17) and (22), the contact constraint equations of each contact pair can be expressed as

Equations (25a) and (25b) are not smooth equations, and the numerical method is not suitable to solve the equations. Thus, it is necessary to deduce the smoothing approximation expression [35] of the contact constraint equations. Equation (25a) can be smoothed and expressed as where is the approximation factor. When approaches , the degree of approximation is extremely high.

The symbol of absolute value is contained in (25b), and it is not convenient to get (25b) smooth via the conventional method. Equation (25b) presents the friction criterion, and it can be illustrated as in Figure 4.

Several functions can be considered to approximate (25b) and the most representative functions can be the following three functions, whose value ranges from −1 to 1:(1)Square root function .(2)Arctangent function .(3)Exponential function .

The expression forms of the square root function and the exponential function are much more complicated than that of the arctangent function. In Section 3.4, the approximated function will be substituted to calculate the dynamic contact force. In order to reduce the computational difficulty, the arctangent function is chosen to approximate (25b). Equation (25b) can be smoothed using the arctangent function [36] and it can be expressed aswhere is the approximation factor. When approaches , the degree of approximation is extremely high. The curves of the approximated step functions, with being , are illustrated in Figure 5.

In general, the allowable embedding distance is rather short and the magnitude of the horizontal coordinate can approach 10⁻⁷. The approximation factor can be reduced to further improve the accuracy. Although the value of is small, and can be calculated within the existing computational capabilities; thus, the operability and accuracy can be ensured.

According to (26) and (27), the approximation expression of (25a) and (25b) can be expressed as

As seen in (18) and (23), the normal relative distance and the tangential relative distance can be expressed with block displacement, while block displacement can be expressed with dynamic contact force , as in (14). Thus, (28a) and (28b) contain only the displacement as an unknown.

Because there are many contact pairs in the block system, every contact pair needs to meet convergence criteria synchronously and they are related to each other. Thus, it is necessary to solve the corresponding (28a) and (28b) for each contact pair simultaneously. Assume that the number of contact pairs is ; contact constraint equations can be expressed as

3.4. Solving Contact Constraint Equations Using Newton’s Method

The convergence rate of the Newton method for solving nonlinear equations is fast and can achieve a square convergence rate; moreover, it has a more stable and more robust calculation convergence [37]. Thus, Newton’s method is adopted to solve the contact constraint equations of the block system . Let be , and the contact constraint equations can be expressed as . Because the contact constraint equations are very complex, the gradient on dynamic contact force cannot be calculated via the analytic method easily. Thus, it is necessary to approximate the gradient using the numerical method. Currently, the difference quotient is widely used to approximate the gradient , and the calculation procedure is as follows.

Let be 0 and specify the initial dynamic contact force . Specify the termination factor .

Calculate the difference quotientwhere controls the accuracy of the calculation and is the unit vector in -dimensional Euclid space.

Solve the following linear equation:

Calculate the error value . Check if the termination condition is verified. If so, stop iterating and set as the optimal solution of the dynamic contact force. If not, conduct the following step.

Set time step as , calculate , and return to step .

After determining the optimal solution of the dynamic contact force , substitute into (14) and the optimal solution of block displacement can be obtained.

4. Boundary Setting

The research area of many dynamic problems is a finite domain, such as the case of impact of the block on the bar in Section 5.1. The boundary of the model can be set as the fixed boundary or the free boundary. Wave reflection occurs when the wave reaches the fixed or free boundary.

The research area of many other dynamic problems is an infinite domain or semi-infinite domain, such as the case of blasting, which is illustrated in Figure 6. The stress wave is expected to propagate freely in the study area. Meanwhile, the actual numerical model used for calculation cannot be infinitely large and the research area can only be a finite domain. If the boundary is not modified, wave reflection occurs when the wave reaches the boundary, which contradicts the real facts [38]. Currently, a viscous boundary is widely used to absorb the energy of the wave propagating to the boundary to avoid wave reflection [39].

The specific method of simulating the viscous boundary is to apply normal viscous force and tangential viscous force at the boundary of the numerical model:where is the density of the block, and are the propagation velocity of the P wave and S wave in the block medium, respectively, and and are the normal and tangential velocity of the blocks at the boundary, respectively.

Based on the theory of wave propagation, and can be determined by solving the following equations:where , , and are Young’s modulus and Poisson’s ratio and density, respectively.

In fact, a viscous boundary is more effective for the internal source problems, such as blasting. It cannot fully simulate earthquakes and other extraneous source problems [40], such as the case of the tunnel under seismic loading in Section 5.2.

Generally, the lateral boundaries are free-field boundaries to simulate earthquakes. Similar to (32), the lateral boundaries are applied using the following equations:where and are the velocities of the block boundary gridpoint in and directions, respectively, and are the velocities of the gridpoint in the free field in and directions, respectively, and and are the force of the gridpoint in the free field in and directions, respectively.

5. Validation

Two cases are chosen to verify the explicit dynamic DDA method considering dynamic contact force proposed in this paper. The case of impact of the block on the bar can verify the correctness of the dynamic contact force of the block under impact dynamic loading. The case of the Xianglushan Tunnel under seismic loading can verify that the new DDA method can be applied in large-scale geotechnical engineering considering complex contact problems.

5.1. Longitudinal Impact of a Block on the Bar

Suppose there is a block with an initial velocity striking towards a bar with the boundary fixed, which is illustrated in Figure 7. The block has the following properties: mass kg and length m. The bar has the following properties: mass kg, density kg/m³, length m, and m. Gravity is not considered in the numerical model. The artificial joints are set in the bar and the bar is divided into 10 subblocks. The value of the artificial joints is large enough to avoid the subblocks being divorced from each other during the process of impact. By setting different initial values of impact velocity and Young’s modulus, the development regularity of dynamic contact force caused by impact of the block on the bar can be studied.

The block and the bar can be seen as a collapsing block system. In the process of impact, the grids of the bar make contact with the moving block, and dynamic contact force is generated in this process. The analytical solution [20] and the numerical solution of the dynamic contact force can be obtained.

In the conventional DDA method, the determination of contact stiffness heavily relies on the engineer’s experience, and the contact stiffness affects the calculation accuracy greatly. Dynamic contact force is introduced in the dynamic DDA method in this paper, and the value of contact stiffness does not have an effect on calculation, but Young’s modulus still affects the calculation accuracy [19]. To study the influence of Young’s modulus on dynamic contact force, the initial impact velocity is taken as m/s and Young’s modulus is taken as , , , and Pa. Figure 8 illustrates the dynamic contact force-time curve with different values of Young’s modulus; the development regularity of the dynamic contact force can also be observed in Figure 8. The dynamic contact force reaches the first peak in a rather short time period after impact and then shows a declining trend. Then the stress wave caused by impacting propagates towards the right fixed boundary in the bar. When the stress wave reaches the right fixed boundary, it can reflect and propagate back towards the left free boundary. The dynamic contact force reaches the second peak when the stress wave reaches the left free boundary. At this time, the time consumed is approximately . Then, the dynamic contact force decreases to 0 within a period time of , where is the length of the bar and is the propagation velocity of the stress wave in the bar. From the aspect of the value of the dynamic contact force, the corresponding value of contact force is larger when the value of Young’s modulus is larger. From the aspect of time period of impact, the corresponding time period of impact is shorter when the value of Young’s modulus is larger, which reveals the rule that the stress wave propagates much faster in the bar with a larger Young’s modulus [41]. Meanwhile, in Figure 8, the dotted line represents the DDA solution and the solid line represents the analytical solution. The analytical solution and the DDA solution of the dynamic contact force are very similar for different values of Young’s modulus, and the dynamic DDA method has high calculation precision.

The initial impact velocity affects the calculation accuracy greatly, and the value of dynamic contact force greatly varies for different initial impact velocities. The initial impact velocity is taken as , and m/s. Figure 9 illustrates a dynamic contact force-time curve for different initial impact velocities. As shown in Figure 9, the development regularity of the dynamic contact force is substantially the same. Two peak values of dynamic contact force are presented in the curve, and the value finally eventually declines to zero. From the aspect of the value of the dynamic contact force, the corresponding value of contact force is larger when the value of initial impact velocity is larger, and the peak value of dynamic contact force is basically proportional to the initial impact velocity. From the aspect of time period of impact, the initial impact velocity has little effect on the impact time, and the development pace of the dynamic contact force for different impact velocities is basically the same. The impact time is only related to the length of the bar and the propagation velocity of the stress wave and has little relation to the initial impact velocity [41]. Meanwhile, the analytical solution and the DDA solution of the dynamic contact force are very similar for different initial impact velocities and the dynamic DDA method has high calculation precision.

According to the analysis of Figures 8 and 9, it can be concluded that the explicit dynamic DDA method considering dynamic contact force proposed in this paper can well describe the dynamic response of contact forces generated by dynamic loading. It can also well reflect the propagation effect of a stress wave in the block.

5.2. Stability Analysis of the Xianglushan Tunnel under Seismic Loading

5.2.1. Project Overview and Numerical Model

Xianglushan Tunnel is in the Yunnan Province of southwest China and is located in a strong earthquake zone with a magnitude of intensity level that reaches VII. The average depth of the tunnel is 700 m, and the maximum depth reaches 900 m. The diameter of the tunnel is 9.8 m, and the whole tunnel spans approximately 62.73 km with several faults. The tunnel is supported by supporting bolts. In the centralized zone of the faults, two groups of widely distributed and fully developed joints exist. The parameters of the joints are N30°~35°E/NW70°~80° and EW/N50°~60°. In this paper, the typical section where the two groups of joints intersect is studied. The DDA model of the Xianglushan Tunnel is illustrated in Figure 10.

The size of the numerical DDA model is m, and the depth of a typical section is 700 m. Because the depth of the tunnel is larger than the size of the tunnel and the number of blocks in the studied area is very large, it is not very convenient to simulate all of the overlying strata above the tunnel. In this case, the numerical DDA model is simplified, being intercepted at 100 m above the tunnel. The size of the tunnel is very small, and the depth is very large; a tunnel with a large depth should be stable compared with a shallow tunnel. The scope of the surrounding rock is more than 10 times larger than the size of the tunnel, and the size of the numerical model can meet the needs of the research [42]. The weight of the overlying strata is applied to the top of the model [43], with the stress of 14.68 MPa. Two groups of joints with a spacing of 3 m are distributed in the model. To facilitate detection, the four typical monitoring points 1, 2, 3, and 4 are set around the tunnel. The mechanical parameters of the surrounding rock and joints are shown in Table 1.

5.2.2. Seismic Loading

Because the underground tunnel is greatly affected by low frequency seismic waves [44, 45] and the KOBE wave contains a rich low frequency band, the KOBE wave is chosen as the input seismic wave, with a duration of 20 s and a peak acceleration of 8.17 m/s². The acceleration-time history curve of the seismic wave needs to be filtered and the baseline corrected. The seismic wave is input from the base of the model. The acceleration-time history curve of the KOBE wave is illustrated in Figure 11.

5.2.3. Calculation Results

In this paper, the displacement-time history of four typical monitoring points 1, 2, 3, and 4 under seismic loading is recorded. The stability of surrounding rock of the tunnel without supporting bolts is studied [46, 47], and the displacement-time history curve of the monitoring points without supporting bolts is illustrated in Figure 12. Because the size of the tunnel is not very large, the displacement-time history curves of monitoring points 1, 2, and 4 are basically coincident and the displacement of monitoring point 3 is slightly larger than the remaining three monitoring points. As observed in Figure 12, without the supporting bolts, the maximum instantaneous displacement of the monitoring points can reach 0.14 m and the permanent displacement after the earthquake can reach 0.014 m.

An implicit cylindric anchor bar element method is adopted for the implementation of the supporting bolts embedded in rock. The embedded supporting bolts are considered to improve the stiffness of the rock in the numerical model. Therefore the stiffness of the supporting bolts can be superimposed onto the stiffness of rock during numerical simulation. This method is detailed in [48].

The displacement-time history curves of the monitoring points with the supporting bolts are illustrated in Figure 13. The displacement-time history curves of monitoring points 1 and 2 are basically coincident. As observed in Figure 13, with the supporting bolts, the maximum instantaneous displacement of the monitoring points can reach 0.056 m and the permanent displacement after the earthquake can reach 0.0068 m. The deformation of the tunnel has been reduced and the stability of the tunnel has been improved significantly compared with the case without supporting bolts.

As can be seen in Figure 13, the displacement of surrounding rock of the Xianglushan Tunnel is small and the anchoring effect on the surrounding rock is obvious. Since the depth of the tunnel is very large, the self-stability of the tunnel can be very good. The two groups of joints existing in the surrounding rock also did not form the blocks which are prone to collapse. Thus, the Xianglushan Tunnel itself can maintain the stability of surrounding rock with the supporting bolts.

In this case, the joints are fully developed and the number of blocks is very large. The contact modes of the blocks are very complex and vary frequently under seismic loading. However, the open-close iteration is avoided in the process of calculation and the contact modes are detected by solving the contact constraint equations of the block system. The contact modes of every time step do not need to be revised repeatedly. As observed in Figures 12 and 13, the results can well reflect the displacement and the deformation of the tunnel under seismic loading, and the calculation result is in accordance with the real case. Thus, we can draw the conclusion that the explicit dynamic DDA method considering dynamic contact force can be applied in the dynamic response analysis of large-scale geotechnical engineering.

6. Discussion

Selection and pairing on contact pairs are often required before contact calculations can be made. When solving the static problems with the conventional static DDA, since the time step is very short, the contact modes of the contact pairs are not modified very frequently. In some continuous time step, the contact modes of the contact pairs can be the same and there is no need to conduct the selection and pairing on contact pairs. Whereas, when solving the dynamic problems with the dynamic DDA, since the dynamic loading acts on the blocks, the contact modes of the contact pairs can be modified very frequently. The repairing of the contact pairs is required at each time step, which increases the computational difficulty. Furthermore, when the number of blocks in the block system is very large and the reciprocating effect of the dynamic loading is very obvious, the selection and pairing of the contact pairs in the early stage can consume much time, which can reduce computational efficiency. Thus, the preselection and the prepairing of the contact pairs under the dynamic loading will be our research focus in the next step.

7. Conclusion

This paper proposed the explicit dynamic DDA method considering dynamic contact force to solve the problem of contact modes detection and dynamic contact force calculation. The calculation process of the explicit dynamic DDA method is deduced and is applied in two numerical cases. Several conclusions can be drawn from this study.

Compared with the conventional DDA method, which adopts block displacement as the unknown for analysis, the dynamic contact force is adopted to solve the problem. The improvement can fully take into account the actual situation instead of taking the dynamic contact force as merely proportional to the embedding distance. This new method can simulate the dynamic contact force between blocks in a more accurate way.

Based on the contact constraint between the block displacement and the dynamic contact force, the contact constraint equations of the block system are deduced, which replaces the role of the open-close iteration in the conventional DDA method. The contact constraint equations can reflect the nonlinear relation of contact constraints and can take into account the interaction between different contact pairs. The contact modes can be obtained once without being revised repeatedly.

The accuracy of the explicit dynamic DDA method is verified by using two numerical cases. In the case of the impact of a block on a bar, taking different values of Young’s modulus and initial impact velocity, the results can well reflect the development regularity of the dynamic contact force. A high degree of agreement exists between the analytical solution and the numerical solution of the dynamic contact force. In the case of stability analysis of the Xianglushan Tunnel under seismic loading, the contact modes of blocks under seismic loading are very complex. The results can reflect the secular deformation behavior of the tunnel under seismic loading. Thus, the explicit dynamic DDA method considering dynamic contact force is still applicable to practical engineering and can be applied in the dynamic response analysis of large-scale geotechnical engineering.

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, Department of Civil Engineering, University of California, Berkeley, Calif, USA, 1988.
X. Pei, R. Huang, and S. Li, “Study of dynamic response of bridge pier shocked by falling rock induced by intensive earthquake,” Chinese Journal of Rock Mechanics and Engineering, vol. 30, no. 2, pp. 3995–4001, 2011.
View at: Google Scholar
J. Ma, X. Zhang, Y. Jiao, and H. Zhang, “Numerical simulation of failure of brazilian disc specimen under dynamic loading using discontinuous deformation analysis method,” Chinese Journal of Rock Mechanics and Engineering, vol. 34, no. 9, pp. 1805–1814, 2015.
View at: Publisher Site | Google Scholar
Y.-J. Ning, J. Yang, and P.-W. Chen, “Numerical simulation of rock blasting in jointed rock mass by DDA method,” Rock & Soil Mechanics, vol. 31, no. 7, pp. 2259–2263, 2010.
View at: Google Scholar
X.-L. Zhang, Y.-Y. Jiao, Q.-S. Liu, and B.-L. Liu, “Numerical study on effect of joints on blasting wave propagation in rock mass,” Rock and Soil Mechanics, vol. 29, no. 3, pp. 717–721, 2008.
View at: Google Scholar
Y. Ning, J. Yang, X. An, and G. Ma, “Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework,” Computers and Geotechnics, vol. 38, no. 1, pp. 40–49, 2011.
View at: Publisher Site | Google Scholar
Y. Ning, J. Yang, G. Ma, and P. Chen, “Modelling rock blasting considering explosion gas penetration using discontinuous deformation analysis,” Rock Mechanics & Rock Engineering, vol. 44, no. 4, pp. 483–490, 2011.
View at: Publisher Site | Google Scholar
J.-H. Wu, “Seismic landslide simulations in discontinuous deformation analysis,” Computers & Geotechnics, vol. 37, no. 5, pp. 594–601, 2010.
View at: Publisher Site | Google Scholar
J.-H. Wu and C.-H. Chen, “Application of DDA to simulate characteristics of the Tsaoling landslide,” Computers and Geotechnics, vol. 38, no. 5, pp. 741–750, 2011.
View at: Publisher Site | Google Scholar
G.-H. Shi, Applications of Discontinuous Deformation Analysis (DDA) to Rock Engineering, Springer, Berlin, Germany, 2007.
X.-J. Kong, J. Liu, and G.-C. Han, “Dynamic failure test and numerical simulation of model concrete-faced rockfill dam,” Chinese Journal of Geotechnical Engineering, vol. 25, no. 1, p. 26, 2003.
View at: Google Scholar
M. Y. Ma, A. D. E. Pan, M. T. Luan et al., “Seismic analysis of stone arch bridges using discontinue deformation analysis,” in Proceedings of the 11th World Conference on Earthquake Engineering, pp. 1–8, June 1996.
View at: Google Scholar
S. Nishiyama, Y. Ohnishi, T. Yanagawa et al., “Study on stability of retaining wall of masonry type by using discontinuous deformation analysis,” in Proceedings of the ISRM International Symposium 3rd ARMS, Y. Ohnishi and K. Aoki, Eds., pp. 1221–1226, Mill Press, Kyoto, Japan, 2004.
View at: Google Scholar
R. Kamai and Y. H. Hatzor, “Numerical analysis of block stone displacements in ancient masonry structures: a new method to estimate historic ground motions,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 32, no. 11, pp. 1321–1340, 2008.
View at: Publisher Site | Google Scholar
C. A. Tang and H. Y. Lu, “The DDD method based on combination of RFPA and DDA. Frontiers of Discontinuous Numerical Methods and Practical Simulations in Engineering and Disaster Prevention,” in Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation (ICADD '11), pp. 105–112, Fukuoka, Japan, August 2013.
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
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 & Geotechnics, vol. 35, no. 3, pp. 346–356, 2008.
View at: Publisher Site | Google Scholar
L.-H. Zhang, T.-Y. Liu, Q.-B. Li, and T. Chen, “Modified dynamic contact force method under reciprocating load,” Engineering Mechanics, vol. 31, no. 7, pp. 8–14, 2014.
View at: Publisher Site | Google Scholar
W. Jiang and H. Zheng, “Influences of artificial parameters in DDA method,” Rock and Soil Mechanics, vol. 28, no. 12, pp. 2603–2606, 2007.
View at: Google Scholar
X. T. Yu, Impact Dynamics, Tsinghua University Press, 2011.
N. Hu, “A solution method for dynamic contact problems,” Computers and Structures, vol. 63, no. 6, pp. 1053–1063, 1997.
View at: Publisher Site | Google Scholar | Zentralblatt MATH
Y. Kanto and G. Yagawa, “A dynamic contact buckling analysis by the penalty finite element method,” International Journal for Numerical Methods in Engineering, vol. 29, no. 4, pp. 755–774, 1990.
View at: Publisher Site | Google Scholar
Y. Kanto and G. Yagawa, “Dynamic contact buckling analysis by the penalty finite element method,” International Journal for Numerical Methods in Engineering, vol. 29, no. 4, pp. 755–774, 1990.
View at: Publisher Site | Google Scholar | Zentralblatt MATH
A. B. Chaudhary and K.-J. Bathe, “A solution method for static and dynamic analysis of three-dimensional contact problems with friction,” Computers and Structures, vol. 24, no. 6, pp. 855–873, 1986.
View at: Publisher Site | Google Scholar
J. C. Simo and T. A. Laursen, “An augmented lagrangian treatment of contact problems involving friction,” Computers & Structures. An International Journal, vol. 42, no. 1, pp. 97–116, 1992.
View at: Publisher Site | Google Scholar | MathSciNet
T. A. Laursen and V. Chawla, “Design of energy conserving algorithms for frictionless dynamic contact problems,” International Journal for Numerical Methods in Engineering, vol. 40, no. 5, pp. 863–886, 1997.
View at: Publisher Site | Google Scholar | MathSciNet
J. Liu, S. Liu, and X. Du, “A method for the analysis of dynamic response of structure containing non-smooth contactable interfaces,” Acta Mechanica Sinic, vol. 15, no. 1, pp. 63–72, 1999.
View at: Publisher Site | Google Scholar
L. Shu, L. Jingbo, S. Jiaguang, and C. Yujian, “A method for analyzing three-dimensional dynamic contact problems in visco-elastic media with kinetic and static friction,” Computers and Structures, vol. 81, no. 24-25, pp. 2383–2394, 2003.
View at: Publisher Site | Google Scholar
Zhou Jing and Ni Hangen, “Response of cracked structures during earthquake,” Earthquake Engineering & Engineering Vibration, vol. 6, no. 2, pp. 52–58, 1986.
View at: Google Scholar
Q. Li, H. Zhou, and G. Lin, “Dynamic fracture analysis of crack at bimaterial interface,” Earthquake Engineering & Engineering Vibration, no. 4, 1990.
View at: Google Scholar
T. JianRen, Dynamic Contact Damping and Its Application in Engineering, Dalian University, 1987.
G. Augusti, “Dynamics of structures: theory and applications to earthquake engineering,” Meccanica, vol. 31, no. 6, pp. 719–720, 1996.
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
H. Zheng and X. Li, “Mixed linear complementarity formulation of discontinuous deformation analysis,” International Journal of Rock Mechanics & Mining Sciences, vol. 75, pp. 23–32, 2015.
View at: Publisher Site | Google Scholar
H.-W. Zhang, S.-Y. He, and X.-S. Li, “Non-interior smoothing algorithm for frictional contact problems,” Applied Mathematics & Mechanics, vol. 25, no. 1, pp. 47–58, 2004.
View at: Publisher Site | Google Scholar
W. Peter and G. Zavarise, Computational Contact Mechanics, Wiley Online Library, 2002.
D. Peter, Newton Methods for Nonlinear Problems, Science Press, 2006.
Y. Zhang, X. Fu, and Q. Sheng, “Modification of the discontinuous deformation analysis method and its application to seismic response analysis of large underground caverns,” Tunnelling and Underground Space Technology, vol. 40, pp. 241–250, 2014.
View at: Publisher Site | Google Scholar
H. Bao, Y. H. Hatzor, and X. Huang, “A new viscous boundary condition in the two-dimensional discontinuous deformation analysis method for wave propagation problems,” Rock Mechanics and Rock Engineering, vol. 45, no. 5, pp. 919–928, 2012.
View at: Publisher Site | Google Scholar
Z. Cui, Q. Sheng, and X. Leng, “Control effect of a large geological discontinuity on the seismic response and stability of underground rock caverns: a case study of the baihetan #1 surge chamber,” Rock Mechanics and Rock Engineering, vol. 49, no. 6, pp. 2099–2114, 2016.
View at: Publisher Site | Google Scholar
J. F. Doyle, Wave Propagation in Structure, Science Press, 2013.
S. G. Qian, X. L. Tang, and L. J. Zhang, “Discussion on calculation limits of 3D subterranean tunnel under dynamic response,” Nonferrous Metals Design, vol. 37, pp. 40–43, 2010.
View at: Google Scholar
Z. Cui, Q. Sheng, and X. Leng, “Effects of a controlling geological discontinuity on the seismic stability of an underground cavern subjected to near-fault ground motions,” Bulletin of Engineering Geology and the Environment, pp. 1–18, 2016.
View at: Publisher Site | Google Scholar
H.-B. Li, X.-D. Ma, W. Shao, X. Xia, H.-C. Dai, and K.-Q. Xiao, “Influence of the earthquake parameters on the displacement of rock cavern,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, pp. 4627–4634, 2005.
View at: Google Scholar
Z. Cui, Q. Sheng, X. Leng, and J. Chen, “Seismic response and stability of underground rock caverns: a case study of Baihetan underground cavern complex,” Journal of the Chinese Institute of Engineers, vol. 39, no. 1, 2016.
View at: Publisher Site | Google Scholar
S. M. Hsiung and G. H. Shi, “Simulation of earthquake effects on underground excavations using Discontinuous Deformation Analysis (DDA),” in Proceedings of the The 38th U.S. Symposium on Rock Mechanics (USRMS '01), ARMA-01-1413, Washington, DC, USA, July 2001.
View at: Google Scholar
Y. H. Hatzor, A. A. Arzi, Y. Zaslavsky, and A. Shapira, “Dynamic stability analysis of jointed rock slopes using the DDA method: king Herod's Palace, Masada, Israel,” International Journal of Rock Mechanics & Mining Sciences, vol. 41, no. 5, pp. 813–832, 2004.
View at: Publisher Site | Google Scholar
Q.-H. Jiang and D.-X. Feng, “Modeling of rockbolts in three-dimensional discontinuous deformation analysis,” Rock and Soil Mechanics, vol. 22, no. 2, pp. 176–178, 2001.
View at: Google Scholar

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

1043

Downloads

2261

Citations