Abstract
The deformation and failure of rock mass is a process of energy dissipation; the damping of DDA is a very important and basic problem. The correctness and effectiveness of DDA rely on the appropriate values of the numeric controlling parameters like time interval, spring stiffness, and assumed maximum displacement ratio , and the contact using the penalty method is the core content of DDA. A mechanical model of two contact blocks loaded with the normal force acting along one side of block boundary is established to study the DDA damping problem, which involves the contact and eliminates the influence of some numeric control parameters (e.g., ). Based on the Newmark method and the theory of DDA, the motion equations of two-block system can be established, and then the relationship of some numeric control parameters and the influence of damping can be obtained. The algorithmic damping increases with the increasing of time interval. Given a very small time interval, the spring stiffness may have no obvious effect on the algorithmic damping. The numerical results reveal that the essence of time interval influencing the open-close iteration is the fact that the algorithmic damping is mainly controlled by time interval.
1. Introduction
The strong earthquakes can destroy rock mass and induce an abundance of landslides, which cause potentially serious threat to both lives and properties. So it is very necessary to study the seismic response of rock mass. Generally, the numerical methods are widely used for dealing with the problem. Furthermore, as the unfavorable geological structure of rock masses, the discontinue-based methods are more suitable for studying the problem.
The discontinuous deformation analysis (DDA), proposed by Shi and Goodman [1, 2], can be used to simulate the discontinuous deformation behavior of jointed rock masses. The correctness and effectiveness of DDA have been verified by MacLaughlin and Doolin [3]. The DDA method is established on the block kinetics; the block dynamic analysis of DDA has been confirmed by many scholars [4–6]. Meanwhile, the DDA method has been widely applied to the seismic response of rock masses. Zhang et al. [7] investigated the seismic response analysis of large underground caverns with the modified DDA method. Zhang et al. [8] studied the run-out analysis of the Daguangbao landslides subjected to near-fault multidirection earthquake forces. Zhang et al. [9] have verified the mobility of earthquake-induced landslide by DDA method. Wu and Chen [10] simulated the kinematic behavior of sliding blocks of rock in the earthquake-induced Tsaoling landslide using seismic discontinuous deformation analysis (DDA). These studies show that DDA method has great potential in the analysis of seismic response of rock masses.
It is generally recognized that the deformation and failure of rock mass is a process of energy dissipation, which should be highlighted in the dynamic problem of rock mass, such as earthquakes. Since the damping research of DDA is the basis of the seismic response, the damping characteristics should be investigated when using DDA to simulate the earthquake-induced dynamic problems of rock mass. In original DDA program, kinetic damping was introduced to discount the initial velocity of each time step. In static analysis the initial velocity is set to be zero for each time step, while the initial velocity for the dynamic analyses has not been discounted. It should be noted that the value of kinetic damping plays a key role in the dynamic analysis of rock mass. Hatzor et al. [11] and Tsesarsky et al. [5] obtained that 2% kinetic damping was the correct number for the dynamic analysis of a single block on an inclined subjected to dynamic loading. However, the significance of the kinetic damping is still not clear, and the value of initial velocity discounted is hard to choose; therefore, some researchers have introduced self-adaptive damping [12] and viscous damping [13, 14] to address these problems. Lin and Xie [15] also analyzed the performance of Newmark time integration with kinetic damping.
Algorithmic damping is caused by the integration scheme of inertia force, which is a controlling factor of the convergence of open-close iteration. Doolin and Sitar [16] studied the time integration of DDA and noted that the algorithmic damping may be important considering the penalty formulation. Bao et al. [17] obtained an equivalent damping ratio from the algorithmic damping in DDA using the cantilever beam model. Overall, the damping researches of DDA are always in progress.
It is well known that the numeric controlling parameters (e.g., assumed maximum displacement ratio , time interval) in DDA are hard to choose; the correctness and effectiveness of DDA depend on the appropriate values of the these numeric control parameters. Moreover, the contact using the penalty method is the core content of DDA, and the open-close iteration convergence is the highlight. Over the years, researches of DDA damping that involve the contact and eliminate the influence of some numeric control parameters (e.g., ) have been rarely reported. In this study, a mechanical model of two contact blocks loaded with the normal force acting along one side of block boundary is established to study the aforementioned problems. Based on the mechanical model, the total equation motion is established easily and intuitively, and the influence of damping is analyzed. The results show that the essence of time interval influencing the open-close iteration is the algorithmic damping mainly controlled by the time interval.
2. Equations of Motion of a Discrete Block System
The DDA method is in essence a numerical method to solve the motion and deformation of multiple rock blocks. Starting from the definitions of displacement functions, the equations of motion of the blocks are established firstly using potential energy minimization. The direct time-integration scheme is adopted to solve the equations of motion.
2.1. Displacement Functions of Deformable Blocks
The DDA method selects the rigid body displacements and the strains of block elements as the basic unknown variables to solve for solutions. Assuming that an arbitrarily shaped polygonal block has uniform stress and strain, an incremental expression for the unknown variables of the block is written aswhere are the translation displacement increments at the point in and directions, is the rotation angle increment of the block around its mass center , and are the increments of normal and shear strains of the block.
Displacement increments of an arbitrary point in the block can be represented aswhere
2.2. Equations of Motion of the Block System Based on the Newmark Method
For a discrete model, the total potential energy of the block system is the sum of all individual blocks’ potential energy. According to the minimum potential principle, the equation of motion of the block system in matrix form [13] can be expressed aswhere is the acceleration vector and , , and are the global mass matrix, stiffness matrix, and load matrix, respectively.
A direct time-integration scheme is used for solving (4), in which the acceleration and the velocity within each time step are taken as the following equations of the Newmark assumption [18]:While denoting , , and , it can be regarded that , , and are the initial velocity vector, the initial acceleration vector, and the displacement increment vector of the block system at the beginning of the current time step, respectively. The kinetic energy dissipation of DDA block system depends on the value of , which is the dynamic controlling parameter. Therefore, is introduced to the Newmark assumption; we haveSubstituting (6) into (4), the motion equation of the block system of the current time step can be derived aswhere
2.3. Algorithmic Damping of the Newmark Method
For analysis of the stability and accuracy of Newmark method, taking , the following equation can be obtained based on (4)-(5): where sampling frequency and undamped frequency of vibration .
The solutions of the characteristic equation about (10) arewhere The real solution must have oscillation in the case of small damping, and the stability must be not infinitely growing; thuswhen is not restricted, (13) should be met, and then , , .
So the conditions of the unconditionally stable Newmark method areIt is thus obvious that [13, 14, 16, 18] both the average acceleration method (, ) and the constant acceleration method (, ) meet the unconditionally stable condition, while the linear acceleration method (, ) and the central difference method (, ) are the conditionally stable method.
In general, the accuracy of Newmark method depends on the time interval, the physical parameters of the system, and the loading condition. Here, we only analyze some widely used unconditionally stable integration schemes, to obtain an idea of the accuracy of these methods evaluated based on (11) . Figure 1 shows the algorithmic damping of the Newmark method versus sampling frequency .
As is shown in Figure 1, of the average acceleration method is unity, which states that this method has no algorithmic damping. By contrast, of the constant acceleration method decreases with the increasing of sampling frequency , which indicates that the constant acceleration method has a rapid attenuation under high frequency interference, while it has very little effect on low frequency. In fact, time interval used in Newmark method is greater than the period of the highest natural frequency of the system; thus, the frequency response is not reliable. So the algorithmic damping is advantageous under certain condition; it may be more suitable for simulating the engineering problems.
3. Two-Block Contact Model and the Submatrices
The DDA method is initially developed for discrete block system. For the fundamental study on the damping of DDA, two-block contact model is established, which involves the contact and aims to eliminate the influence of some numeric control parameters (e.g., ).
3.1. Contact Constraint Conditions of DDA
The DDA method can simulate large displacements of blocks and allows blocks to slide, rotate, and open. A key of the DDA method is contact problem; within each time step, no penetration and no tension between the block interfaces must be satisfied. For 2D DDA method, the contact types include vertex-to-vertex, vertex-to-edge, and edge-to-edge contacts. The edge-to-edge contact can usually be converted to two vertex-to-edge contacts in DDA original program. In this paper, the contact of two blocks is locked, and a normal spring and a shear spring are applied on the contact interface of the two blocks.
3.2. Two-Block Contact Model
The model of two contact blocks loaded with the normal force applied along one side of block boundary as shown in Figure 2 is established to analyze the damping problem quantitatively. In the model, block B is totally fixed, and block A is subjected to the uniform distributed load on the upper boundary. The length of block A is and the gravity as well as the contact friction is not considered.
Based on the Newmark method and the theory of DDA, the motion equations of two-block system can be established, and then the relationship of some parameters and the influence of damping can be obtained.
3.3. Calculations of the Submatrices
As block B is totally fixed, the deformation matrix of block B should be taken as . The deformation matrix of block A denotes . It is clear that , and the initial deformation matrix of block A is . So the global mass matrix, stiffness matrix, and load matrix of two-block system are given as , , and .
In the performance of DDA calculation, several submatrices, such as the inertia force submatrix, elastic submatrix, spring stiffness submatrix, load submatrix, and initial stress submatrix, should be deduced and assembled.
3.3.1. Inertia Force Submatrix
The mass matrix of block A is where is the unit mass.where is the area of block . , , and , where is the integral of over block, is the integral of over block, is the integral of over block, is the integral of over block, and is the integral of over block.
For block A, , , , , , , , , and . Therefore, (15) can be rewritten asThe corresponding inertia force submatrix can be obtained as
3.3.2. Elastic Submatrix
The strain energy of block A iswhere is the expanded elastic matrix.
The corresponding elastic submatrix can be obtained by minimizing the strain energy where and are the elastic modulus and the poisson ratio, respectively.
3.3.3. Spring Stiffness Submatrix
The stiffness submatrix of normal spring [2] iswhere , , and As shown in Figure 2, , , and are , , and , respectively. Consequently is simplified asSimilarly, for the contact point , , , and are , , and , respectively. Thus of block is .
The stiffness submatrix [2] of shear spring iswhereSimilarly, , .
As a result, the spring stiffness submatrix of block can be acquired as
3.3.4. Load Submatrix
As is shown in Figure 2, the uniform load is equal to . The corresponding load submatrix [2] of block can be presented as
3.3.5. Initial Stress Submatrix
The corresponding initial stress submatrix is where , , and are the stress components of th time step of block .
4. Damping Analysis of DDA
Based on the two-block contact model, the influencing factors of damping can be analyzed quantitatively. By assembling the global stiffness matrix and the load matrix, the total equation motion can be obtained.
4.1. Formula of the Total Equation Motion
Solving (8) based on two-block model, is formed by (18), (21), and (27), is formed by (19) and (28)–(30). As , the following three iteration equations can be obtained:where, , and are the deformation increments at the th time step of block , and are the initial velocity and initial acceleration at the th time step of block , and are the accumulated displacement at points and , and and are the stress components at the th time step.
After the values , , and are determined, using the initial conditions, , , , and , , , and can be obtained by solving (31). Substituting into (2), the physical equations of elastic dynamics, (6) and (7), respectively, , , , and can be obtained.
4.2. Algorithmic Damping of DDA
4.2.1. Influencing Factors of Algorithmic Damping
For studying the algorithmic damping and eliminating the influence of kinetic damping, the kinetic damping is not considered by denoting . The acceleration within each time step of DDA is assumed to be constant (, ).
Considering that , , , and are the physical and mechanical parameters, we only study the influences of time interval and the spring stiffness . In order to perform case studies, assumptions are made that m, kg/m3, GPa, , and MPa. The curve of displacement of the middle point of block versus time step under the different time interval (0.01 sec, 0.002 sec, 0.001 sec) is shown in Figure 3. The corresponding contact force curve is presented in Figure 4. The curve of displacement of the middle point of block versus time step under the different spring stiffness (1, 10, 100) is shown in Figure 5. The corresponding contact force curve is shown in Figure 6. The results show that contact force approximates the real solution 5 MPa successively and keeps stable at last, which satisfies the equilibrium of force system.
As shown in Figures 3 and 4, the smaller the time interval is, the more powerful the displacement and contact force oscillate. With the increasing of time interval, the displacement and contact force can converge to the equilibrium of position rapidly, which is actually caused by the algorithmic damping. For example, if sec, the displacement and contact force reach their convergence state after five and four time steps, respectively. But if sec, the time consumption for convergence of displacement and contact force is thirty-four time steps and thirty-six time steps, respectively, much more than those of sec.
Clearly, from Figures 5 and 6, as the value of spring stiffness decreases, the oscillation becomes more obvious. However, when spring stiffness increases to a certain degree, it causes no difference. For example, if , the convergence of both displacement and contact force costs forty-four time steps. By contrast, if , the convergence of displacement and contact force costs twenty and twenty-one time steps, respectively. In addition, the difference in convergence of and is not obvious.
Furthermore, when time interval is 0.00001 sec, the relationship between displacement of the middle point and time steps tends to be vibration curves with no attenuation, and spring stiffness only influences the amplitude and frequency of vibration curves, which can be seen in Figure 7.
4.2.2. Analysis of Algorithmic Damping
The deformation and failure of rocks is a process of energy dissipation; the damping research of DDA is a fundamental problem. Algorithm damping can be caused by DDA using Newmark constant acceleration integral scheme. If average acceleration method is used in DDA simulation, there will be no algorithmic damping, which can be seen in Figure 8. Significantly, the algorithm damping is useful under certain conditions, it can guarantee the open-close convergence, and for the quasi-static problem, it can quickly converge to the correct solution.
From Section 4.2.1, it can be seen that time interval is a critical factor of the algorithmic damping. It has been shown that, with increasing the time interval, the displacement and contact force can converge to the equilibrium of position rapidly. However, time interval is an artificial parameter because the values have significant effect on the accuracy of the results. For instance, if the time interval is too large, it may cause the loss of diagonal dominance, interfering with the accuracy. If the time interval is too small, the algorithmic damping may be neglected; thus the contact displacement may vibrate indefinitely, which seriously affects the convergence of open-close iteration.
Moreover, the smaller the value of spring stiffness is, the more obvious the oscillation is, and when spring stiffness is large enough, it causes no difference. These phenomena illustrate that the algorithmic damping decreases with the decreasing of spring stiffness under certain conditions. Additionally, when the time interval is too small, the spring stiffness may have no obvious effect on the algorithmic damping and only influence the amplitude and frequency of vibration curve, which explains that time interval is more sensitive to the algorithmic damping.
4.3. Kinetic Damping of DDA
For a better understanding of the kinetic damping, the acceleration within each time step of DDA is assumed to be average acceleration where , ; in that case the algorithmic damping of DDA can be neglected. The physical and mechanical parameters are the same as Section 4.2.1. The curve of displacement of the middle point of block versus time steps under the different (1.00, 0.98, 0.95, 0.90) is shown in Figure 9.
Some scholars have researched the kinetic damping, and most of them made their efforts to compare the simulation results with physical-experimental results, which aimed to obtain the discounted values of the initial velocity of each time step. However, because of the complexity of the energy dissipation, the accurate discounted value is still difficult to determine. Here, we investigate the relative trend of kinetic damping under the different values of dynamic control parameter, where the corresponding discounted values are 0%, 2%, 5%, and 1%. Taking sec and , curves of the displacement against time steps can be obtained as given in Figure 9, from which the amplitude of vibration curve is found to attenuate fast with the increase of the discounted value.
5. Conclusions
The analysis of the seismic response is of great significance in rock engineering. The DDA method is very suitable for the problem, which has been confirmed by many scholars. However, the deformation and failure of rocks is a process of energy dissipation; the damping research of DDA is a fundamental problem. In this paper, we perform quantitative analysis of the influencing factors of damping based on the Newmark method and the mechanical model of two contact blocks. Through the analysis of damping, the following understandings are obtained.
() The time interval is a controlling factor of the algorithmic damping. The algorithmic damping increases with the increasing of time interval; the displacement and contact force can converge to the equilibrium of position rapidly with larger time interval.
() The spring stiffness has a certain effect on the algorithmic damping. The algorithmic damping increases as spring stiffness increases, but if spring stiffness is large enough, it will cause no difference.
() Given a small enough time interval, the algorithmic damping may be neglected, and the spring stiffness may have no obvious effect on the algorithmic damping; it can only influence the amplitude and frequency of vibration curve. It is clear that the time interval is more sensitive to algorithmic damping. The essence of time interval influencing the open-close iteration is the algorithmic damping mainly controlled by the time interval.
() The kinetic damping of DDA is only dependent on the discounted values of the initial velocity of each time step. The amplitude of vibration curve attenuates fast with the increase of the discounted value, which illustrates that the appropriate discounted value is critical to the energy dissipation.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank Dr. G. H. Shi for his valuable advice in understanding of DDA method. Financial support for this work, provided by the National Natural Science Foundation of China (no. 51679173) and Natural Science Foundation of Hubei Province (no. 2016CFA083), is gratefully acknowledged.