# Explicit Dynamic DDA Method considering Dynamic Contact Force

Authors: Jian Zhao; Ming Xiao; Juntao Chen; Dongdong Li
Journal: Shock and Vibration
 (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.


---

## Body

## 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 blocki contacts the adjacent blocks; the contact effect between the block i and the adjacent blocks can be regarded as the point loading, which is illustrated in Figure 1.Figure 1
Contact model of blocks.The contact force applied on blocki as point loading can be classified into two categories. The first category is active contact force Ri1, as illustrated in Figure 1, and it is applied to the contact point of block i itself. The second category is passive contact force Ri2, as illustrated in Figure 1, and it is applied to the contact line of block i itself. From the principle of interacting force, the passive contact force Ri2 applied to block i and the active contact force Rj2 applied to block j are a pair of interacting forces; thus, the relationship between the two interacting forces can be expressed as(1)Ri2=-Rj2.The active contact forceRi1 can be decomposed into RXi1 in the direction of x and RYi1 in the direction of y, and the contact potential energy generated by active contact force applied to block i can be expressed as (2)Πp=-RXi1u1+RYi1v1=-u1,v1RXi1RYi1=-DiTTi1x,yTRXi1RYi1,where Di is the displacement of block i and Ti1(x,y) is the displacement function at the contact point of block i. Taking the derivative of -Πp on displacement Di, a 6×1 submatrix of point loading can be obtained and expressed as (3)Ti1x,yTRXi1RYi1=Ti1x,yTRi1⟶Ri.Meanwhile, the point loading submatrix corresponding to passive contact forceRi2 can be expressed as (4)Ti2x,yTRXi2RYi2=Ti2x,yTRi2⟶Ri.To unify the context, all of the equations are deduced with active contact force. With reference to (1), (4) can be expressed as (5)Ti2x,yTRi2=-Ti2x,yTRj2⟶Ri.After the contact submatrix of blocki is calculated, the remaining submatrices, including elastic submatrix, initial stress submatrix, loading submatrix, and displacement-constraint submatrix, except the inertia force submatrix of block i, are calculated and integrated into the stiffness matrix Ki and the loading matrix Fi. With reference to (4) and (5), the simultaneous equations of block i can be expressed as (6)KiUi=Fi+Ri.Assume that the number of blocks of the whole block system isNb. A similar deduction can be made for every block, and the simultaneous equations can be expressed as(7)KiUi=Fi+Rii=1,…,Nb.Integrating (7) of all blocks, the overall simultaneous equations of the block system can be expressed as(8)KU=F+R.All 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 submatrixMU¨ into (8), the overall simultaneous equations of the block system considering inertia force can be expressed as (9)MU¨+KU=F+R,where M is the overall mass matrix of the block system after being integrated and the mass submatrix of each single block Mi can be expressed as (10)Mi=ρ∬ΩiTiTTidxdyi=1⋯Nb,where ρ is the density of a block and Ti 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(11)U¨t=U¨t-Δt=2Δt2Ut-2ΔtU˙t-Δt.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 as(12)U¨t=1Δt2Ut-Δt-2Ut+Ut+Δt.Substituting (12) into (9), the discretization of the overall simultaneous equations in the time domain can be expressed as (13)M1Δt2Ut-Δt-2Ut+Ut+Δt+KUt=Ft+Rt.The displacement at time step t+Δt can be obtained using (13) and can be expressed as(14)Ut+Δt=U-t+Δt+Δt2M-1Rt,U-t+Δt=2Ut-Ut-Δt+Δt2M-1Ft-KUt.It can be seen from (14) that the displacement U(t+Δt) at time step t+Δt is determined by the displacement at time steps t and t-Δt, which is the distinctive feature of the explicit algorithm [32]. From (14), the displacement U(t+Δt) at time step t+Δt consists of two parts [18]: the predictable displacement U-(t+Δt) without considering dynamic contact force and the amended displacement Δt2M-1R(t) caused by dynamic contact force. The predictable displacement U-(t+Δt) is only related to the displacement and the external load at time steps t and t-Δt and is a known quantity. Thus, the dynamic contact force R(t) needs to be solved to obtain the displacement at time step t+Δt.


## 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 blocki contacts the adjacent blocks; the contact effect between the block i and the adjacent blocks can be regarded as the point loading, which is illustrated in Figure 1.Figure 1
Contact model of blocks.The contact force applied on blocki as point loading can be classified into two categories. The first category is active contact force Ri1, as illustrated in Figure 1, and it is applied to the contact point of block i itself. The second category is passive contact force Ri2, as illustrated in Figure 1, and it is applied to the contact line of block i itself. From the principle of interacting force, the passive contact force Ri2 applied to block i and the active contact force Rj2 applied to block j are a pair of interacting forces; thus, the relationship between the two interacting forces can be expressed as(1)Ri2=-Rj2.The active contact forceRi1 can be decomposed into RXi1 in the direction of x and RYi1 in the direction of y, and the contact potential energy generated by active contact force applied to block i can be expressed as (2)Πp=-RXi1u1+RYi1v1=-u1,v1RXi1RYi1=-DiTTi1x,yTRXi1RYi1,where Di is the displacement of block i and Ti1(x,y) is the displacement function at the contact point of block i. Taking the derivative of -Πp on displacement Di, a 6×1 submatrix of point loading can be obtained and expressed as (3)Ti1x,yTRXi1RYi1=Ti1x,yTRi1⟶Ri.Meanwhile, the point loading submatrix corresponding to passive contact forceRi2 can be expressed as (4)Ti2x,yTRXi2RYi2=Ti2x,yTRi2⟶Ri.To unify the context, all of the equations are deduced with active contact force. With reference to (1), (4) can be expressed as (5)Ti2x,yTRi2=-Ti2x,yTRj2⟶Ri.After the contact submatrix of blocki is calculated, the remaining submatrices, including elastic submatrix, initial stress submatrix, loading submatrix, and displacement-constraint submatrix, except the inertia force submatrix of block i, are calculated and integrated into the stiffness matrix Ki and the loading matrix Fi. With reference to (4) and (5), the simultaneous equations of block i can be expressed as (6)KiUi=Fi+Ri.Assume that the number of blocks of the whole block system isNb. A similar deduction can be made for every block, and the simultaneous equations can be expressed as(7)KiUi=Fi+Rii=1,…,Nb.Integrating (7) of all blocks, the overall simultaneous equations of the block system can be expressed as(8)KU=F+R.All 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 submatrixMU¨ into (8), the overall simultaneous equations of the block system considering inertia force can be expressed as (9)MU¨+KU=F+R,where M is the overall mass matrix of the block system after being integrated and the mass submatrix of each single block Mi can be expressed as (10)Mi=ρ∬ΩiTiTTidxdyi=1⋯Nb,where ρ is the density of a block and Ti 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(11)U¨t=U¨t-Δt=2Δt2Ut-2ΔtU˙t-Δt.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 as(12)U¨t=1Δt2Ut-Δt-2Ut+Ut+Δt.Substituting (12) into (9), the discretization of the overall simultaneous equations in the time domain can be expressed as (13)M1Δt2Ut-Δt-2Ut+Ut+Δt+KUt=Ft+Rt.The displacement at time step t+Δt can be obtained using (13) and can be expressed as(14)Ut+Δt=U-t+Δt+Δt2M-1Rt,U-t+Δt=2Ut-Ut-Δt+Δt2M-1Ft-KUt.It can be seen from (14) that the displacement U(t+Δt) at time step t+Δt is determined by the displacement at time steps t and t-Δt, which is the distinctive feature of the explicit algorithm [32]. From (14), the displacement U(t+Δt) at time step t+Δt consists of two parts [18]: the predictable displacement U-(t+Δt) without considering dynamic contact force and the amended displacement Δt2M-1R(t) caused by dynamic contact force. The predictable displacement U-(t+Δt) is only related to the displacement and the external load at time steps t and t-Δt and is a known quantity. Thus, the dynamic contact force R(t) needs to be solved to obtain the displacement at time step t+Δt.


## 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 Figure2.Figure 2
Model of normal contact constraint.When the contact mode is embedding, (15) is met and expressed as(15)dn=0,Rn>0.When the contact mode is open, (16) is met and expressed as (16)dn>0,Rn=0,where dn is the normal relative distance and Rn is the normal contact force. According to (15) and (16), (17) can be expressed as follows:(17)mindn,Rn=0.Normal relative distancedn can be expressed as (18)dn=δn+n→TjDj-TiDi,where δn is the initial spacing of the contact pair at the beginning of the time step, Dj and Di are the displacements of the two contacted blocks, respectively, Tj and Ti are the displacement matrices of the two nearest points of the contact pair at the beginning of the time step, and n→ 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(19)Rn=n→R.


### 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 Figure3.Figure 3
Model of tangential contact constraint.When the contact mode is sliding, (20) is met and expressed as (20)dτ=0,μRn-Rτ>0.When the contact mode is fixed, (21) is met and expressed as (21)dτ>0,μRn-Rτ=0,where dτ is the tangential relative distance, Rn is the normal contact force, Rτ is the tangential contact force, and μ is the friction coefficient. According to (20) and (21), (22) can be expressed as(22)min⁡dτ,μRn-Rτ=0.Tangential relative distancedτ can be expressed as (23)dτ=τ→TjDj-TiDi,where Dj and Di are the displacements of the two contacted blocks, respectively, Tj and Ti 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 forceRτ in the local coordinate system and the contact force R in the global coordinate system. The coordinate transform can be expressed as(24)Rτ=τ→R.


### 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 (25a)mindn,Rn=0,(25b)mindτ,μRn-Rτ=0.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 (26)-ε1ln⁡exp⁡-dnε1+exp⁡-Rnε1=0,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.Figure 4
Frictional constraint condition.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 functionG1Δ=Δ/Δ2+ε.(2)
Arctangent functionG2Δ=(2/π)arctanΔ/ε.(3)
Exponential functionG3Δ=1-2/(eΔ/ε+1).The expression forms of the square root function and the exponential function are much more complicated than that of the arctangent function. In Section3.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(27)Rτ-2πμRnarctan⁡dτε2=0,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.Figure 5
Approximated function curve of the step function.In general, the allowable embedding distancedτ 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, dτ 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(28a)-ε1ln⁡exp⁡-dnε1+exp⁡-Rnε1=0,(28b)Rτ-2πμRnarctan⁡dτε2=0.As seen in (18) and (23), the normal relative distance dn and the tangential relative distance dτ can be expressed with block displacement, while block displacement D can be expressed with dynamic contact force R, as in (14). Thus, (28a) and (28b) contain only the displacement D 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 Nc; contact constraint equations HRn,Rτ can be expressed as (29)HRn,Rτ=-ε1lnexp-dniε1+exp-Rniε1Rτi-2πμRniarctandτiε2i=Nci=1=0.


### 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 HRn,Rτ. Let Rn,Rτ be R, and the contact constraint equations can be expressed as HR. Because the contact constraint equations are very complex, the gradient on dynamic contact force R 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 VR, and the calculation procedure is as follows.( 1 ) Let k be 0 and specify the initial dynamic contact force R0. Specify the termination factor ε>0.( 2 ) Calculate the difference quotient(30)VkRk=HRk+hke-HRkhk,where hk controls the accuracy of the calculation and e is the unit vector in n-dimensional Euclid space.( 3 ) Solve the following linear equation:(31)VkRkΔRk+HkRk=0.( 4 ) Calculate the error value gRk=HRkTHRk. Check if the termination condition gRk≤ε′ is verified. If so, stop iterating and set Rk as the optimal solution of the dynamic contact force. If not, conduct the following step.( 5 ) Set time step k as k+1, calculate Rk+1=Rk+ΔRk, and return to step (2).( 6 ) After determining the optimal solution of the dynamic contact force R, substitute R into (14) and the optimal solution of block displacement D can be obtained.


## 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 Figure2.Figure 2
Model of normal contact constraint.When the contact mode is embedding, (15) is met and expressed as(15)dn=0,Rn>0.When the contact mode is open, (16) is met and expressed as (16)dn>0,Rn=0,where dn is the normal relative distance and Rn is the normal contact force. According to (15) and (16), (17) can be expressed as follows:(17)mindn,Rn=0.Normal relative distancedn can be expressed as (18)dn=δn+n→TjDj-TiDi,where δn is the initial spacing of the contact pair at the beginning of the time step, Dj and Di are the displacements of the two contacted blocks, respectively, Tj and Ti are the displacement matrices of the two nearest points of the contact pair at the beginning of the time step, and n→ 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(19)Rn=n→R.


## 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 Figure3.Figure 3
Model of tangential contact constraint.When the contact mode is sliding, (20) is met and expressed as (20)dτ=0,μRn-Rτ>0.When the contact mode is fixed, (21) is met and expressed as (21)dτ>0,μRn-Rτ=0,where dτ is the tangential relative distance, Rn is the normal contact force, Rτ is the tangential contact force, and μ is the friction coefficient. According to (20) and (21), (22) can be expressed as(22)min⁡dτ,μRn-Rτ=0.Tangential relative distancedτ can be expressed as (23)dτ=τ→TjDj-TiDi,where Dj and Di are the displacements of the two contacted blocks, respectively, Tj and Ti 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 forceRτ in the local coordinate system and the contact force R in the global coordinate system. The coordinate transform can be expressed as(24)Rτ=τ→R.


## 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 (25a)mindn,Rn=0,(25b)mindτ,μRn-Rτ=0.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 (26)-ε1ln⁡exp⁡-dnε1+exp⁡-Rnε1=0,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.Figure 4
Frictional constraint condition.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 functionG1Δ=Δ/Δ2+ε.(2)
Arctangent functionG2Δ=(2/π)arctanΔ/ε.(3)
Exponential functionG3Δ=1-2/(eΔ/ε+1).The expression forms of the square root function and the exponential function are much more complicated than that of the arctangent function. In Section3.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(27)Rτ-2πμRnarctan⁡dτε2=0,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.Figure 5
Approximated function curve of the step function.In general, the allowable embedding distancedτ 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, dτ 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(28a)-ε1ln⁡exp⁡-dnε1+exp⁡-Rnε1=0,(28b)Rτ-2πμRnarctan⁡dτε2=0.As seen in (18) and (23), the normal relative distance dn and the tangential relative distance dτ can be expressed with block displacement, while block displacement D can be expressed with dynamic contact force R, as in (14). Thus, (28a) and (28b) contain only the displacement D 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 Nc; contact constraint equations HRn,Rτ can be expressed as (29)HRn,Rτ=-ε1lnexp-dniε1+exp-Rniε1Rτi-2πμRniarctandτiε2i=Nci=1=0.


## 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 HRn,Rτ. Let Rn,Rτ be R, and the contact constraint equations can be expressed as HR. Because the contact constraint equations are very complex, the gradient on dynamic contact force R 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 VR, and the calculation procedure is as follows.( 1 ) Let k be 0 and specify the initial dynamic contact force R0. Specify the termination factor ε>0.( 2 ) Calculate the difference quotient(30)VkRk=HRk+hke-HRkhk,where hk controls the accuracy of the calculation and e is the unit vector in n-dimensional Euclid space.( 3 ) Solve the following linear equation:(31)VkRkΔRk+HkRk=0.( 4 ) Calculate the error value gRk=HRkTHRk. Check if the termination condition gRk≤ε′ is verified. If so, stop iterating and set Rk as the optimal solution of the dynamic contact force. If not, conduct the following step.( 5 ) Set time step k as k+1, calculate Rk+1=Rk+ΔRk, and return to step (2).( 6 ) After determining the optimal solution of the dynamic contact force R, substitute R into (14) and the optimal solution of block displacement D 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 Section5.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 Figure6. 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].Figure 6
DDA model under blasting loading.The specific method of simulating the viscous boundary is to apply normal viscous forcefn and tangential viscous force fs at the boundary of the numerical model:(32)fn=ρcpvn,fs=ρcsvs,where ρ is the density of the block, cp and cs are the propagation velocity of the P wave and S wave in the block medium, respectively, and vn and vs are the normal and tangential velocity of the blocks at the boundary, respectively.Based on the theory of wave propagation,cp and cs can be determined by solving the following equations:(33)cp=E1-νρ1+ν1-2ν,cs=E2ρ1+ν,where E, ν, 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:(34)fx=-ρcpvxm-vxffS+fxff,fy=-ρcsvym-vyffS+fyff,where vxm and vym are the velocities of the block boundary gridpoint in x and y directions, respectively, vxff and vyff are the velocities of the gridpoint in the free field in x and y directions, respectively, and fxff and fyff are the force of the gridpoint in the free field in x and y 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 velocityv0 striking towards a bar with the boundary fixed, which is illustrated in Figure 7. The block has the following properties: mass M1=1.0 kg and length L1=0.075 m. The bar has the following properties: mass M2=1.5 kg, density ρ=2700 kg/m3, length L2=0.5 m, and L3=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.Figure 7
DDA model of impact of a block on the bar.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 v0=12.5 m/s and Young’s modulus is taken as E=0.56×107, 1.0×107, 2.25×107, and 9.0×107 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 2L2/c. Then, the dynamic contact force decreases to 0 within a period time of L2/c, where L2 is the length of the bar and c 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.Figure 8
Dynamic contact force-time curve with different values of Young’s modulus.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 asv0=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.Figure 9
Dynamic contact force-time curve for different initial impact velocities.According to the analysis of Figures8 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 Figure10.Figure 10
DDA model of Xianglushan Tunnel under seismic loading.The size of the numerical DDA model is200×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.Table 1
Mechanical parameters of rock and joints.
Materials Density (kg/m3) Young’s modulus (GPa) Poisson’s ratio Friction angle (°) Cohesion (MPa) Tensile strength (MPa) Rock 2850 5 0.29 — — — Joints — — — 35 0.4 0


#### 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/s2. 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.Figure 11
Acceleration-time curve of the KOBE wave.


#### 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.Figure 12
Displacement-time curve of monitoring points without supporting bolts.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 Figure13. 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.Figure 13
Displacement-time curve of monitoring points with supporting bolts.As can be seen in Figure13, 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 Figures12 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.


## 5.1. Longitudinal Impact of a Block on the Bar

Suppose there is a block with an initial velocityv0 striking towards a bar with the boundary fixed, which is illustrated in Figure 7. The block has the following properties: mass M1=1.0 kg and length L1=0.075 m. The bar has the following properties: mass M2=1.5 kg, density ρ=2700 kg/m3, length L2=0.5 m, and L3=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.Figure 7
DDA model of impact of a block on the bar.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 v0=12.5 m/s and Young’s modulus is taken as E=0.56×107, 1.0×107, 2.25×107, and 9.0×107 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 2L2/c. Then, the dynamic contact force decreases to 0 within a period time of L2/c, where L2 is the length of the bar and c 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.Figure 8
Dynamic contact force-time curve with different values of Young’s modulus.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 asv0=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.Figure 9
Dynamic contact force-time curve for different initial impact velocities.According to the analysis of Figures8 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 Figure10.Figure 10
DDA model of Xianglushan Tunnel under seismic loading.The size of the numerical DDA model is200×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.Table 1
Mechanical parameters of rock and joints.
Materials Density (kg/m3) Young’s modulus (GPa) Poisson’s ratio Friction angle (°) Cohesion (MPa) Tensile strength (MPa) Rock 2850 5 0.29 — — — Joints — — — 35 0.4 0


### 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/s2. 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.Figure 11
Acceleration-time curve of the KOBE wave.


### 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.Figure 12
Displacement-time curve of monitoring points without supporting bolts.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 Figure13. 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.Figure 13
Displacement-time curve of monitoring points with supporting bolts.As can be seen in Figure13, 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 Figures12 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.


## 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 Figure10.Figure 10
DDA model of Xianglushan Tunnel under seismic loading.The size of the numerical DDA model is200×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.Table 1
Mechanical parameters of rock and joints.
Materials Density (kg/m3) Young’s modulus (GPa) Poisson’s ratio Friction angle (°) Cohesion (MPa) Tensile strength (MPa) Rock 2850 5 0.29 — — — Joints — — — 35 0.4 0


## 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/s2. 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.Figure 11
Acceleration-time curve of the KOBE wave.


## 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.Figure 12
Displacement-time curve of monitoring points without supporting bolts.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 Figure13. 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.Figure 13
Displacement-time curve of monitoring points with supporting bolts.As can be seen in Figure13, 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 Figures12 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.( 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.



---

*Source: 7431245-2016-12-28.xml*