Explicit Dynamic DDA Method considering Dynamic Contact Force

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.


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-erhshan 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 openclose 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.

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.

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  1 , as illustrated in Figure 1, and it is applied to the contact point of block  itself.The second category is passive contact force  2 , 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  2 applied to block  and the active contact force  2 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  1 can be decomposed into  1 in the direction of  and  1 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 [ 1 (, )] is the displacement function at the contact point of block .Taking the derivative of −Π  on displacement [  ], a 6 × 1 submatrix of point loading can be obtained and expressed as Meanwhile, the point loading submatrix corresponding to passive contact force  2 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 displacementconstraint 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 as Integrating ( 7) of all blocks, the overall simultaneous equations of the block system can be expressed as All of the contact force in ( 8) is taken as active contact force, and inertia force is not taken into consideration.

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 openclose 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 as Substituting ( 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 (Δ) 2 [] −1 {()} 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  + Δ.

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.

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 as When 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

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) 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  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  1 is the approximation factor.When  1 approaches 0 + , 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  1 (Δ) = Δ/ √ Δ 2 + .
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 as where  2 is the approximation factor.When  2 approaches 0 + , the degree of approximation is extremely high.The curves of the approximated step functions, with  2 being 6.0 × 10 −8 , are illustrated in Figure 5.
In general, the allowable embedding distance   is rather short and the magnitude of the horizontal coordinate Δ can approach 10 −7 .The approximation factor  2 can be reduced to further improve the accuracy.Although the value of  2 is small,   and  2 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

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.(1) Let  be 0 and specify the initial dynamic contact force { (0) }.Specify the termination factor  > 0.
(2) Calculate the difference quotient where ℎ  controls the accuracy of the calculation and {} is the unit vector in -dimensional Euclid space.
(3) Solve the following linear equation: (4) 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.
(6) After determining the optimal solution of the dynamic contact force {}, substitute {} into (14) and the optimal solution of block displacement {} can be obtained.

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 p widely used to absorb the energy of the wave propagating to the boundary to avoid wave reflection [39].

Viscous boundary
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 V  and V  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 V   and V   are the velocities of the block boundary gridpoint in  and  directions, respectively, V   and V   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.

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.

Longitudinal Impact of a
Block on the Bar.Suppose there is a block with an initial velocity V 0 striking towards a bar with the boundary fixed, which is illustrated in Figure 7.The block has the following properties: mass  1 = 1.0 kg and length  1 = 0.075 m.The bar has the following properties: mass  2 = 1.5 kg, density  = 2700 kg/m 3 , length  2 = 0.5 m, and  3 = 0.05 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 V 0 = 12.5 m/s and Young's modulus is taken as  = 0.56 × 10 7 , 1.0 × 10 7 , 2.25 × 10 7 , and 9.0 × 10 7 Pa. Figure 8 illustrates the dynamic contact forcetime 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 2 2 /.Then, the dynamic  contact force decreases to 0 within a period time of  2 /, where  2 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 V 0 = 12.5, 25.0, 37.5, and 50.0 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.

Stability Analysis of the Xianglushan Tunnel under
Seismic Loading

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 200 × 300 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.

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 2 .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.

Calculation Results.
In this paper, the displacementtime 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 largescale geotechnical engineering.

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.

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.
(1) 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.
(2) 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.
(3) 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.

Figure 1 :
Figure 1: Contact model of blocks.

Figure 2 :
Figure 2: Model of normal contact constraint.

Figure 3 :
Figure 3: Model of tangential contact constraint.

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

1 Figure 7 :
Figure 7: DDA model of impact of a block on the bar.

Figure 8 :
Figure 8: Dynamic contact force-time curve with different values of Young's modulus.

Figure 9 :
Figure 9: Dynamic contact force-time curve for different initial impact velocities.

Figure 12 :Figure 13 :
Figure 12: Displacement-time curve of monitoring points without supporting bolts.

Table 1 :
Mechanical parameters of rock and joints.