Investigation of 2D Seismic DDA Method for Numerical Simulation of Shaking Table Test of Rock Mass Engineering

Abstract

Since the basic theory of the discontinue deformation analysis (DDA) method was proposed, the DDA open source has gone through a long development process. At present, different kinds of programs have been widely applied in rock mass engineering such as slope, dam, and tunnel. This paper introduces the solution principle of DDA motion equations in detail, as well as the development status of the 2D open-source program. Numerical simulation of shaking table test of rock mass engineering using 2D DDA program is highlighted, and investigations of seismic wave pre-processing and seismic input method are carried out. First, based on the Newmark integration scheme, the integration algorithms of synthetic or measured seismic wave time history, correction function of seismic wave, and DDA simulation are unified. Then, three seismic input methods are implanted in the DDA program, and the applicability of various seismic input methods is discussed. On this basis, using the improved seismic 2D DDA program, a shaking table test of typical rock mass engineering is simulated. Through the comparison between the theoretical/test data and simulation results, the reliability of the improved DDA program in seismic response analysis is verified; the large mass method and the large stiffness method are more suitable for rigid foundation, such as shaking table test; the propagation of the seismic wave presents a significant amplification effect due to the reflection, refraction, and diffraction in the tunnel. The research results provide DDA theory and an open-source program for analyzing the seismic response of rock mass engineering.

1. Introduction

In the past century, a series of major earthquakes have occurred all over the world. The rock mass damage caused by earthquakes has attracted great attention in the engineering field [1]. Therefore, the dynamic response analysis of rock mass engineering under seismic action is an important subject to be broken through. Usually, the seismic dynamic response of rock mass engineering mostly adopts the finite element method, finite difference method, finite volume method, boundary element method, etc. [2,3,4,5,6,7,8]. These numerical methods concentrate on the elastic-plastic analysis. There are few open-source programs that can reproduce the failure process. To solve this problem, many scholars have carried out theoretical research and program development. At present, the methods based on discontinuous medium theory, such as discrete element method and discontinuous deformation analysis (DDA) method, are being developed rapidly. The DDA method is a new numerical simulation technique for geotechnical stability analysis, and it was originally proposed by Shi G.H. in his doctoral thesis [9,10,11]. DDA is a method parallel to the finite element method. The difference is that the DDA method can not only reflect the displacement and deformation of a block, but also allow the block to have motion modes such as sliding, opening, and rotation.

Since it was proposed, the DDA method has been deeply explored. Since 1996, the international conference on analysis of discontinuous deformation has been held regularly, and 15 sessions have been held so far. On the correctness and effectiveness of DDA method, Ma, Maclaughlin and Yagoda-Biran et al. [12,13,14] carried out a comprehensive review of relevant articles. However, the original DDA method also has limitations in dealing with many problems, such as high-precision stress field in a rock block, stratified rock mass, and so on. In order to solve these problems, many scholars have carried out in-depth research, and the theoretical research can be roughly summarized as parameter selection, block contact, accurate block stress, dynamic calculation, multiscale analysis, and 3D simulation [15,16]. Meanwhile, the DDA method has also been widely used in practical engineering, including slopes, dams, tunnels, and so on [17,18].

The DDA method adopts the Newmark integration scheme to simulate the discontinuous motion process of the rock block system in real time-steps. This approach makes it a potential method for analyzing the seismic response of rock mass engineering. In slope analysis, Wu, Zhang, Fu, Huang, and Wang et al. [19,20,21,22,23] simulated the failure process of starting, moving, and stacking of Tsaoling landslide, Daguangbao landslide, Ganjiazhai landslide, Niumiangou landslide, and Donghekou landslide under seismic action. In terms of dam analysis, Shyu, Kong, and Kaidi et al. [24,25,26] evaluated the seismic stability of a typical arch dam and rockfill dam. Relatively speaking, there are few DDA results of the seismic response of tunnels [27,28]. In addition, the DDA method is also applied to the seismic response analysis of landfill, bridge, retaining wall, and masonry structures [29,30,31,32,33].

Although DDA has been applied to the seismic response analysis of a number of rock mass engineering examples, the verification of field monitoring data is still insufficient. Compared with the scarcity of field monitoring data, the measured data based on shaking table test is easier to obtain. Akao and Sasaki et al. [34,35] simulated the seismic response of single block and three block systems using DDA, which was in good agreement with the test results. Shimaoka, Irie, and Feng et al. [36,37,38] reproduced the shaking table test of a rock slope by DDA simulation, and the comparison between the test data and simulation result proved that the DDA method is reasonable to simulate seismic response of rock mass engineering. It should be noted that at present, the 2D DDA open-source program is mainly used in analyzing the seismic response of the rock mass engineering, while the theory and program of the 3D DDA method are still developing [39,40,41]. Therefore, the 2D DDA method is used in this paper. Based on the existing 2D DDA applied to numerical simulation of shaking table test, this paper introduces the development status of the open-source program of 2D DDA in detail, the seismic wave pre-processing and seismic input method are investigated to extend the 2D DDA program; using the seismic DDA method, a shaking table test of tunnel engineering is simulated, and the propagation law of seismic wave in rock mass is analyzed.

2. A Brief Description of 2D DDA Method and Its Open-Source Program

2.1. Motion Equations Based on Newmark Integration Scheme

The motion equations of DDA method are used to describe the response of block system under load, and its basis is the Lagrange function:

$$L\left(\left\{D\right\},\left\{\dot{D}\right\}\right)=T\left(\left\{\dot{D}\right\}\right)-V\left(\left\{D\right\}\right)$$

(1)

where {Ḋ} and {D} are the velocity vector and displacement vector, respectively; T() and V() are functions of kinetic energy and potential energy, respectively. In the block system simulated by DDA, the potential energy is mainly composed of elastic stress, initial stress, point load, volume load, and point displacement of a block, as well as the contacts between blocks.

The Lagrange function is time integrated, and then the motion equations are obtained by the Hamilton principle:

$$\frac{d}{dt}\frac{\partial T}{\partial \left\{\dot{D}\right\}}-\frac{\partial V}{\partial \left\{D\right\}}=0$$

(2)

Expanding the above equation, the dynamic ordinary differential equations of the block system can be obtained,

$$\left[M\right]\left\{\ddot{D}\right\}+\left[C\right]\left\{\dot{D}\right\}+\left[K\right]\left\{D\right\}=\left\{F\right\}$$

(3)

where $\left\{\ddot{D}\right\}$ is the acceleration vector; [M], [C], [K], and {F} are the mass matrix, damping matrix, stiffness matrix, and load vector of the block system, respectively.

The solution methods of motion equations mainly include mode superposition method and direct integration method. The solution object of DDA is a discrete block system. Once the sliding and separation between block elements occurs, its vibration mode changes constantly, so it is difficult to use the mode superposition method. The direct integration method does not need to change the equation form before solving the motion equation. It directly carries out numerical integration step-by-step, and is more suitable for solving the differential equations of DDA.

The solution process of the direct integration method is based on two concepts. First, it is assumed that the time solution domain [0, t] of the motion equations is equally divided into N time-steps, and in each time-step, the motion equations should be satisfied. Second, in a time-step, the displacement, velocity, and acceleration meet an approximate functional expression. Therefore, according to Equation (3), the motion equations of the block system in time-step n + 1 are:

$$\left[M\right]{\left\{\ddot{D}\right\}}_{n+1}+\left[C\right]{\left\{\dot{D}\right\}}_{n+1}+\left[K\right]{\left\{D\right\}}_{n+1}={\left\{F\right\}}_{n+1}$$

(4)

According to different expressions of the approximate function, the direct integration method is divided into the Houbolt integration scheme, Wilson integration scheme, and Newmark integration scheme. The Newmark integration scheme is commonly used by engineers [42,43]. The relationships among displacement, velocity, and acceleration described by the Newmark integration scheme are as follows:

$${\left\{D\right\}}_{n+1}={\left\{D\right\}}_n+\Delta t{\left\{\dot{D}\right\}}_n+\frac{\Delta t^2}{2}\left[\left(1-2\beta \right)\left\{\ddot{D}\right\}{}_n+2\beta \left\{\ddot{D}\right\}{}_{n+1}\right]$$

(5)

$${\left\{\dot{D}\right\}}_{n+1}={\left\{\dot{D}\right\}}_n+\Delta t\left[\left(1-\gamma \right)\left\{\ddot{D}\right\}{}_n+\gamma \left\{\ddot{D}\right\}{}_{n+1}\right]$$

(6)

where Δt is a time-step; β and γ are integration parameters. In the open-source program released by Dr. Shi G.H., γ = 1 and β = 1/2, which is the expression of the constant acceleration integration method. Then, substituting Equations (5) and (6) into Equation (4), the equilibrium equations of the block system are obtained:

$$\left[\widehat{K}\right]{\left\{D\right\}}_{n+1}={\left\{\widehat{F}\right\}}_{{}_{n+1}}$$

(7)

where,

$$\begin{array}{c}\left\{\widehat{K}\right\}=\left[K\right]+\frac{1}{\beta \Delta t^2}\left[M\right]+\frac{\gamma }{\beta \Delta t}\left[C\right]\\ {\left\{\widehat{F}\right\}}_{{}_{n+1}}={\left\{F\right\}}_{n+1}+\left[M\right]\left[\frac{1}{\beta \Delta t^2}{\left\{D\right\}}_n+\frac{1}{\beta \Delta t}{\left\{\dot{D}\right\}}_n+\left(\frac{1}{2\beta }-1\right){\left\{\ddot{D}\right\}}_n\right]+\left[C\right]\left[\frac{\gamma }{\beta \Delta t}{\left\{D\right\}}_n+\left(\frac{\gamma }{\beta }-1\right){\left\{\dot{D}\right\}}_n+\left(\frac{\gamma }{2\beta }-1\right)\Delta t{\left\{\ddot{D}\right\}}_n\right]\end{array}$$

In the DDA method, the displacement of each deformable discrete block is a solving target. For a 2D problem, six mechanical parameters with obvious physical significance are selected as the incremental displacement of a block element in one time-step. Taking block-i as an object,

$$\Delta D_i={\left\{u_0,v_0,r_0,{\epsilon }_x,{\epsilon }_y,{\gamma }_{xy}\right\}}^T$$

(8)

where {ΔD_i} is the solution displacement of the block; (u₀, v₀) is the block centroid for the rigid body translation; r₀ is the rotation angle around the block centroid; and ε_x, ε_y, and γ_xy are the normal and shear strains of the block deformation.

For the first order assumption, the incremental displacements of any point (x, y) of block-i in the current time-step can be obtained:

$$\left\{\begin{array}{c}\Delta u_i\\ \Delta v_i\end{array}\right\}=\left[T_i\right]\Delta D_i$$

(9)

where Δu_i and Δv_i are components of the incremental displacements in the x direction and the y direction, respectively; $\left[T_i\right]=\left[\begin{array}{cccccc}1& 0& -\left(y-y_0\right)& x-x_0& 0& \left(y-y_0\right)/2\\ 0& 1& x-x_0& 0& -\left(y-y_0\right)& \left(x-x_0\right)/2\end{array}\right]$ is the shape function matrix, and (x₀, y₀) is the centroid of block-i.

Supposing there are n blocks in the defined block system, based on the Equation (7), the equilibrium equations of are written as:

$$\left[\begin{array}{cccc}{{\displaystyle K}}_{11}& {{\displaystyle K}}_{12}& \cdots & {{\displaystyle K}}_{1n}\\ {{\displaystyle K}}_{21}& {{\displaystyle K}}_{22}& \cdots & {{\displaystyle K}}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {{\displaystyle K}}_{n1}& {{\displaystyle K}}_{n2}& \cdots & {{\displaystyle K}}_{nn}\end{array}\right]\left\{\begin{array}{c}{{\displaystyle \Delta D}}_1\\ {{\displaystyle \Delta D}}_2\\ ⋮\\ {{\displaystyle \Delta D}}_n\end{array}\right\}=\left\{\begin{array}{c}{{\displaystyle F}}_1\\ {{\displaystyle F}}_2\\ ⋮\\ {{\displaystyle F}}_n\end{array}\right\}$$

(10)

2.2. Open-Source Program of 2D DDA

Since the basic theory of 2D DDA was put forward, the open-source program has gone through a long development process. Among all developers, Dr. Shi G.H., the initiator of the DDA method, made the greatest contribution. In 1986, the first DDA code was developed on PC with Basic language and was written by Dr. Shi G.H. In 1989, Dr. Shi G.H. used NDP C compiler to complete the first version of DDA code with C language. Since then, he has modified the DDA code many times with C language. In order to improve the generality of DDA, many scholars have also developed or improved the open-source program, as shown in

Table 1. Existing open-source program of 2D DDA.

Code Version	Language	Developer	Remarks
DDA’1986	BASIC	Shi G.H.	The first DDA code
DDA’1989	C/PC	Shi G.H.	The first C DDA code
DDA’1992	C/UNIX	Shi G.H.	The first X-windows DDA code
DDA’1994	C/UNIX	Shi G.H.	Cohesion, tensile strength and auto stiffness
DDA’1995	C/UNIX	Shi G.H.	Correction of shear lock and rotation problem
DDA/W’1995	C/Win	MacLaughlin and Doolin [44]	The first DDA/Win DDA code
DDA/WT’1995	C/SUN	Ohnishi et al. [45]	Advanced, new functions, only on SUN
DDA/AJ’1994	FORTRAN	Ke and Goodman [46]	Artificial joint
DDA/SB’1995	FORTRAN	Lin et al. [47]	Sub-block, breakable block
DDAML’2001	XML	Doolin and Sitar [48]	Processing input files with XML
HDDA’2008	C, C++/Win	He [49]	Input CAD; enrich post-processing function
CDDA’2009	C/Win	Zheng and Jiang [50]	Complementary theory
DDA-SC’2010	C/Win	Khan [51]	Deformation joint model
DDA/Zhao’2012	C++/Win	Zhao et al. [52]	GUI based on shell concept
DDARF’2012	C, C++/Win	Jiao et al. [53]	Boundary or internal cracking algorithm
easyDDA’2013	C, C++/Win	Cheng et al. [54]	An easy-to-use GUI software

.

The open-source program of 2D DDA used in this paper was provided by Dr. Shi G.H. in 2012. This DDA code was developed based on Windows and programmed with C language. The program consists of four modules, including DDA Lines (DL), DDA Cut (DC), DDA Forward (DF), and DDA Graph (DG). DL module uses the Monte Carlo method to generate a joint network. The main function of the DC module is to form a block system based on computational geometry and combinatorial topology theory. DF is the core module of the DDA program. It solves the displacement and stress of a block under external load by inputting geometric model, physical and mechanical parameters, seismic load, and other data. Its analysis process is shown in Figure 1. DG is a post-processing module. It can draw the displacement and stress vector of a block and reproduce the failure process of the block system.

3. Seismic Wave Pre-Processing Based on Newmark Integration Scheme for DDA Program

3.1. Pre-processing Method

As shown in Figure 1, before the DDA program is executed, input data of seismic load are needed. Generally, the seismic wave observed by seismic stations or synthesized is acceleration data. However, the input data of seismic wave of DDA simulation can be acceleration, velocity, or displacement. Therefore, it is necessary to integrate the acceleration once to obtain the velocity, and the displacement is obtained by secondary integration. During these integration processes, there may be a drift phenomenon, which needs baseline correction [55]. Most correction functions are continuous functions, and there are integrals in the process of correcting discrete seismic wave data.

As shown in Figure 2, synthetic or measured seismic wave integration and baseline correction of seismic wave are pre-processed to provide reliable input data for DDA simulation. It is necessary to ensure the unity of the integration scheme of the pre-processing and the DDA simulation. In the DDA method, the Newmark integration scheme is used. Therefore, this integration scheme is also required for seismic wave integration and baseline correction.

It is assumed that the acceleration data observed by seismic stations or synthesized are known, (a₀, a₁, a₂,…, a_m−1, a_m). The time-step of the data is Δt, the total duration is T, the initial seismic velocity is V₀ = 0, and the initial displacement d₀ = 0. Using the Newmark integration scheme, the velocity can be obtained by integrating the acceleration once:

$$v_{n+1}=v_n+\Delta t\left[\left(1-\gamma \right)a{}_n+\gamma a{}_{n+1}\right] (n=0, 1, \dots , \mathrm{m}-1)$$

(11)

The displacement can be obtained by integrating the acceleration twice:

$$d_{n+1}=d_n+\Delta t{v}_n+\frac{\Delta t^2}{2}\left[\left(1-2\beta \right)a{}_n+2\beta a{}_{n+1}\right] (n=0, 1, \dots , \mathrm{m}-1)$$

(12)

At the end of seismic action, the acceleration a_m may be non-zero, or v_m and d_m are not equal to zero after integration, which means that there is a drift phenomenon, and the baseline correction is required. It is assumed that the time domain of baseline correction is [0, t], and the correction functions of acceleration, velocity, and displacement are expressed as I_a(t), I_v(t), and I_d(t), respectively. For the correction functions, the following two conditions must be met [56,57]:

(a)

The initial values of acceleration, velocity, and displacement must be zero.

(b)

The termination values of acceleration, velocity, and displacement must be equal to −a_m, −v_m, and −d_m, respectively, which are the opposite of the termination value of seismic action.

Generally, there are six unknowns in the above two conditions, and six linear equations need to be solved. However, if the appropriate correction function is selected to automatically meet the first condition (a), the unknowns can be reduced to 3. In this paper, the cubic polynomial is selected as the correction function of acceleration, and its expression is:

$$I_a\left(t\right)=at+b{t}^2+c{t}^3$$

(13)

where a, b, and c are the unknown coefficients.

On the premise that the initial value is zero, the first and quadratic integral expressions of I_a(t) are:

$$I_v\left(t\right)=\frac{1}{2}a{t}^2+\frac{1}{3}b{t}^3+\frac{1}{4}c{t}^4$$

(14)

$$I_d\left(t\right)=\frac{1}{6}a{t}^3+\frac{1}{12}b{t}^4+\frac{1}{20}c{t}^5$$

(15)

Equations (13)–(15) automatically satisfy the first condition (a). According the second condition (b):

$$\begin{array}{c}aT+b{T}^2+c{T}^3+a_m=0\\ \frac{1}{2}a{T}^2+\frac{1}{3}b{T}^3+\frac{1}{4}c{T}^4+v_m=0\\ \frac{1}{6}a{T}^3+\frac{1}{12}b{T}^4+\frac{1}{20}c{T}^5+d_m=0\end{array}$$

(16)

The unique values of a, b, and c can be obtained by solving Equation (16), and the correction functions of acceleration, velocity, and displacement can be obtained by substituting them into Equations (13)–(15). In the process of obtaining Equations (14) and (15) by integrating Equation (13), the accurate integration of continuous function is adopted. In fact, due to the seismic data being composed of a number of discrete points and the DDA simulation also being in discrete time-step, the correction function should be integrated in the form of discrete points. On this basis, a baseline correction method based on Newmark integration scheme is proposed.

For the cubic polynomial of the correction function, the seismic acceleration at each time-step can be expressed as:

$$I_a\left(\mathrm{n}\right)=an\Delta t+b{\left(n\Delta t\right)}^2+c{\left(n\Delta t\right)}^3$$

(17)

The baseline correction method based on Newmark integration scheme includes the following four steps.

Step 1: The velocity correction function is obtained by integrating the acceleration correction function:

$$I_v\left(\mathrm{n}+1\right)={\displaystyle \underset{k=0}{\overset{k=n+1}{\int }}{I}_a\left(\mathrm{k}\right)d\tau =}{\displaystyle \underset{k=0}{\overset{k=n}{\int }}{I}_a\left(k\right)d\tau +}{\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}{I}_a\left(k\right)d\tau =}{I}_v\left(n\right)+{\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}{I}_a\left(k\right)d\tau }$$

(18)

Using the relationship between the velocity and the acceleration described by Newmark integration scheme,

$$\begin{array}{cc}\hfill {\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}{I}_a\left(k\right)d\tau }\approx & \Delta t\left[\left(1-\gamma \right)I_a\left(n\right)+\gamma I_a\left(n+1\right)\right]=\left(n+\gamma \right)\Delta t^{2}a\hfill \\ \hfill +& \left(n^2+2\gamma n+\gamma \right)\Delta t^{3}b+\left(n^3+3\gamma n^2+3\gamma n+\gamma \right)\Delta t^{4}c\hfill \end{array}$$

(19)

Substituting Equation (19) into Equation (18), the velocity correction function is written as:

$$\begin{array}{l}{I}_v\left(\mathrm{n}\right)=I_v\left(0\right)+{\displaystyle \sum _{i=0}^{i=n-1}\left({\displaystyle \underset{k=i}{\overset{k=i+1}{\int }}{I}_a\left(\mathrm{k}\right)d\tau }\right)}\approx I_v\left(0\right)+\\ \begin{array}{ccccccc}& & & & & & {\displaystyle \sum _{i=0}^{i=n-1}\left[\left(i+\gamma \right)\Delta t^{2}a+\left(i^2+2\gamma i+\gamma \right)\Delta t^{3}b+\left(i^3+3\gamma i^2+3\gamma i+\gamma \right)\Delta t^{4}c\right]}\end{array}\end{array}$$

(20)

Step 2: The displacement correction function is obtained from the correction function of acceleration and velocity:

$$\begin{array}{l}{I}_{\mathrm{d}}\left(\mathrm{n}+1\right)={\displaystyle \underset{k=0}{\overset{k=n+1}{\int }}{I}_v\left(\mathrm{k}\right)d\tau =}{\displaystyle \underset{k=0}{\overset{k=n}{\int }}{I}_v\left(k\right)d\tau +}{\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}{I}_v\left(k\right)d\tau =}{I}_d\left(n\right)+{\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}{I}_v\left(k\right)d\tau }\\ \begin{array}{cccccccccccc}& & & & & & & & & & & \begin{array}{ccc}\begin{array}{ccc}& & \end{array}& & =I_d\left(n\right)+I_v\left(n\right)\Delta t+{\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}\left[{\displaystyle \underset{i=k-1}{\overset{i=k}{\int }}{I}_a\left(i\right)d\tau }\right]d\tau }\end{array}\end{array}\end{array}$$

(21)

Using the relationship between the displacement and the acceleration described by Newmark integration scheme,

$$\begin{array}{l}{\displaystyle \underset{k=n}{\overset{k=n+1}{\int }}\left[{\displaystyle \underset{i=k-1}{\overset{i=k}{\int }}{I}_a\left(i\right)d\tau }\right]d\tau }\approx \left(\frac{1}{2}-\beta \right)\Delta t^2{I}_a\left(n\right)+\beta \Delta t^2{I}_a\left(n+1\right)\\ =\frac{\Delta t^2}{2}\left[\left(1-2\beta \right)I_a\left(n\right)+2\beta I_a\left(n+1\right)\right]\\ =\left[\left(n+2\beta \right)\frac{\Delta t^{3}a}{2}+\left(n^2+4\beta n+2\beta \right)\frac{\Delta t^{4}b}{2}+\left(n^3+6\beta n^2+6\beta n+2\beta \right)\frac{\Delta t^{5}c}{2}\right]\end{array}$$

(22)

Substituting Equation (22) into Equation (21), the displacement correction function is written as:

$$\begin{array}{l}{I}_d\left(\mathrm{n}\right)\approx I_d\left(0\right)+{\displaystyle \sum _{i=0}^{i=n-1}\left(I_v\left(0\right)\Delta t+{\displaystyle \sum _{j=0}^{j=i-1}\left[\left(j+\gamma \right)\Delta t^{3}a+\left(j^2+2\gamma j+\gamma \right)\Delta t^{4}b+\left(j^3+3\gamma j^2+3\gamma j+\gamma \right)\Delta t^{5}c\right]}\right)}\\ +{\displaystyle \sum _{i=0}^{i=n-1}\left[\left(i+2\beta \right)\frac{\Delta t^{3}a}{2}+\left(i^2+4\beta i+2\beta \right)\frac{\Delta t^{4}b}{2}+\left(i^3+6\beta i^2+6\beta i+2\beta \right)\frac{\Delta t^{5}c}{2}\right]}\end{array}$$

(23)

Step 3: Combined with Equations (17), (20) and (23), the values of parameters a, b, and c can be calculated by solving Equation (16), and the acceleration, velocity, and displacement correction functions are obtained.

Step 4: Using the correction functions, the acceleration, velocity, and displacement at each time-step are corrected.

3.2. Example Verification

Taking a synthetic seismic wave as an example, its acceleration time history is shown as the initial wave in Figure 3a. The duration of the seismic action is 40 s and the time-step is 0.02 s [28]. Using the same Newmark integration scheme as in the DDA program, where γ = 1 and β = 1/2, the acceleration is integrated into velocity and displacement, which are shown as the initial wave in Figure 3b,c. It is not difficult to find that there is a drift phenomenon in the integration of seismic data. Using the proposed method for baseline correction, the correction waves in Figure 3 are obtained, which show that the drift phenomenon has been well eliminated. Specifically,

Table 2. Initial and corrected values at different time-steps.

Time/s	Acceleration/(m/s2)		Velocity/(m/s)		Displacement/m
	Initial	Correction	Initial	Correction	Initial	Correction
0	0	0	0	0	0	0
5	−0.114560	−0.112080	0.081587	0.089092	0.025958	0.039160
10	0.432663	0.435038	−0.088690	−0.068190	−0.227670	−0.144980
15	−0.01504	−0.014340	−0.139280	−0.110660	−0.175540	0.032970
20	−0.03556	−0.037090	−0.105020	−0.078490	−0.144660	0.206531
25	0.184509	0.181186	0.027189	0.041156	−0.198760	0.258175
30	0.448024	0.444370	−0.086800	−0.091100	−0.207480	0.275434
35	−0.054380	−0.05590	0.019649	0.001206	−0.391630	0.030864
40	−0.004090	0	0.013603	0	−0.330420	0

lists the initial and corrected values at different time-steps. At the end of seismic action (40 s), the values of initial acceleration, velocity, and displacement are −0.004090 m/s², 0.013603 m/s, and −0.330420 m, respectively. After correcting based on the Newmark integration scheme, all values are 0. Moreover, with the progress of integration, the correction magnitude of the initial value becomes larger and larger, which means the correction magnitude meets this law: displacement > velocity > acceleration.

4. Seismic Input Method in the DDA Program

4.1. Multi-Blocks Newmark Method

In the original DDA program written by Dr. Shi G.H., the seismic loading is applied to all blocks in the form of inertial force, and this treatment extends the classical Newmark method to the block system, so we name it the multi-blocks Newmark method.

By the coordinate transformation and the linear interpolation of the synthetic or measured seismic wave, the seismic accelerations for the DDA simulation in each time-step are obtained, and the inertial force acted on the unit area of block-i is:

$$\left\{\begin{array}{c}{f}_x\\ f_y\end{array}\right\}=\rho \left\{\begin{array}{c}{a}_x\\ a_y\end{array}\right\}$$

(24)

where f_x and f_y are the components of the inertial force in the x direction and the y direction, respectively, ρ is the block medium density, and a_x and a_y are the seismic accelerations in the x direction and the y direction, respectively.

Then, the potential energy of the inertia force is:

$${\Pi }_q=-{\displaystyle \iint \left\{\begin{array}{cc}\Delta u_i& \Delta v_i\end{array}\right\}\left\{\begin{array}{c}{f}_x\\ f_y\end{array}\right\}dxdy}=-\Delta D_i{}^T{\displaystyle \iint {\left[T_i\right]}^{T}dxdy}\rho \left\{\begin{array}{c}{a}_x\\ a_y\end{array}\right\}$$

(25)

where ${\displaystyle \iint {\left[T_i\right]}^{T}dxdy}=\left[\begin{array}{cc}{S}_i& 0\\ 0& S_i\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]$ and S_i is the block area.

Using the principle of minimum potential energy, the contributions of the seismic loading to the equilibrium equations of DDA are obtained:

$$\rho S_i\left\{\begin{array}{c}{a}_x\\ a_y\\ 0\\ 0\\ 0\\ 0\end{array}\right\}\to \left\{F_i\right\}$$

(26)

4.2. Large Mass Method

The large mass method is another acceleration input method. In this method, the bottom block of the DDA model is special and has a large mass but no gravity [58]; the seismic loading is applied to the bottom block and is timely passed to the block system.

The governing equations of the block system are rewritten as:

$$\left[M\right]\left\{\ddot{D}\right\}+\left[C\right]\left\{\dot{D}\right\}+\left[K\right]\left\{D\right\}=\left[M_{bb}\right]\left\{{\ddot{D}}_g\right\}$$

(27)

where [M_bb] is the mass matrix of the bottom block, $\left\{{\ddot{D}}_g\right\}$ is the seismic acceleration vector.

It can be decomposed into another form:

$$\left[\begin{array}{cc}{M}_{ss}& 0\\ 0& M_{bb}\end{array}\right]\left\{\begin{array}{c}\ddot{D}{}_s\\ {\ddot{D}}_b\end{array}\right\}+\left[\begin{array}{cc}{C}_{ss}& 0\\ 0& C_{bb}\end{array}\right]\left\{\begin{array}{c}\dot{D}{}_s\\ {\dot{D}}_b\end{array}\right\}+\left[\begin{array}{cc}{K}_{ss}& K_{sb}\\ K_{sb}{}^T& K_{bb}\end{array}\right]\left\{\begin{array}{c}D{}_s\\ D_b\end{array}\right\}=\left\{\begin{array}{c}0\\ M_{bb}{\ddot{D}}_g\end{array}\right\}$$

(28)

where [M_ss] and [C_ss] are the mass matrix and the damping matrix, respectively, of the block system except for the bottom block; [C_bb] is the damping matrix of the bottom block; the subscript s and b represent the block system except for the bottom block and the bottom block, respectively.

Expanding the second row of Equation (28):

$$\left[M_{bb}\right]\left\{{\ddot{D}}_b\right\}+\left[C_{bb}\right]\left\{{\dot{D}}_b\right\}+{\left[K_{sb}\right]}^T\left\{D{}_s\right\}+\left[K_{bb}\right]D_b=\left[M_{bb}\right]\left\{{\ddot{D}}_g\right\}$$

(29)

Both sides of Equation (29) are pre-multiplied by [M_bb]⁻¹ and it becomes:

$$\left\{{\ddot{D}}_b\right\}+{\left[M_{bb}\right]}^{-1}\left[C_{bb}\right]\left\{{\dot{D}}_b\right\}+{\left[M_{bb}\right]}^{-1}{\left[K_{sb}\right]}^T\left\{D{}_s\right\}+{\left[M_{bb}\right]}^{-1}\left[K_{bb}\right]D_b=\left\{{\ddot{D}}_g\right\}$$

(30)

As a result of the large mass of the bottom block, diagonal elements of [M_bb]⁻¹ tend to zero, if the damping of the bottom block is independent of the mass, Equation (29) can be approximately written as:

$${\left\{\ddot{D}\right\}}_b\approx {\left\{\ddot{D}\right\}}_g$$

(31)

Thus, the large mass method is a kind of approximate solution, and it inputs a seismic acceleration by setting a large number in the mass matrix. Equation (26) can be used to calculate the contributions of the seismic loading of the bottom block to the equilibrium equations of DDA.

4.3. Large Stiffness Method

The large stiffness method is a displacement input method. In this method, rigid springs are forced to the bottom block of the DDA model; the seismic loading is applied to the rigid springs in the form of displacement and is timely passed to the block system.

The governing equations of the block system are rewritten as:

$$\left[M\right]\left\{\ddot{D}\right\}+\left[C\right]\left\{\dot{D}\right\}+\left[K\right]\left\{D\right\}=\left[K_{bb}\right]\left\{D_g\right\}$$

(32)

where [K_bb] is the stiffness matrix of the bottom block, $\left\{D_g\right\}$ is the seismic displacement vector.

It can be decomposed into another form:

$$\left[\begin{array}{cc}{M}_{ss}& 0\\ 0& M_{bb}\end{array}\right]\left\{\begin{array}{c}\ddot{D}{}_s\\ {\ddot{D}}_b\end{array}\right\}+\left[\begin{array}{cc}{C}_{ss}& 0\\ 0& C_{bb}\end{array}\right]\left\{\begin{array}{c}\dot{D}{}_s\\ {\dot{D}}_b\end{array}\right\}+\left[\begin{array}{cc}{K}_{ss}& K_{sb}\\ K_{sb}{}^T& K_{bb}\end{array}\right]\left\{\begin{array}{c}D{}_s\\ D_b\end{array}\right\}=\left\{\begin{array}{c}0\\ K_{bb}{D}_g\end{array}\right\}$$

(33)

Expanding the second row of Equation (33):

$$\left[M_{bb}\right]\left\{{\ddot{D}}_b\right\}+\left[C_{bb}\right]\left\{{\dot{D}}_b\right\}+{\left[K_{sb}\right]}^T\left\{D{}_s\right\}+\left[K_{bb}\right]D_b=\left[K_{bb}\right]\left\{D_g\right\}$$

(34)

Both sides of Equation (34) are pre-multiplied by [K_bb]⁻¹ and it becomes:

$${\left[K_{bb}\right]}^{-1}\left[M_{bb}\right]\left\{{\ddot{D}}_b\right\}+{\left[K_{bb}\right]}^{-1}\left[C_{bb}\right]\left\{{\dot{D}}_b\right\}+{\left[K_{bb}\right]}^{-1}{\left[K_{sb}\right]}^T\left\{D{}_s\right\}+D_b=\left\{D_g\right\}$$

(35)

As a result of the rigid springs of the bottom block, diagonal elements of [K_bb]⁻¹ tend to zero, if the damping of the bottom block is independent of the stiffness, Equation (35) can be approximately written as:

$$\left\{D_b\right\}\approx \left\{D_g\right\}$$

(36)

Thus, the large stiffness method is also a kind of approximate solution. It is very convenient to implement this method in the DDA program by the displacement load boundary, and the pseudo code is shown in Figure 4.

4.4. Example Verification

Based on the original 2D DDA code, using the basic theories of the above three seismic input methods, the seismic DDA program is extended. At the same time, using the proposed baseline correction method in Section 3, the input waves required by the seismic DDA program are pre-processed. In order to test the correctness of three different seismic input methods in 2D DDA program, a verification example was designed.

The geometry of the established DDA model is shown in Figure 5a, where the size of the fixed block A is 16 m × 4 m and the size of block B is 4 m × 4 m. The mechanical parameters of the rock include an elastic Young’s modulus of 10 GPa, a Poisson’s ratio of 0.25, a density of 2500 kg/m³, and a gravity acceleration of 10 m/s². The mechanical parameters of the contact face include a cohesion strength of zero, a tensile strength of zero, a friction angle of 30°, and a spring coefficient of 2.5 GPa.

As shown in Figure 6, the input dynamic loading is horizontal and its acceleration time history is a sine wave of a period a(t) = 10 sin(πt), where 0 ≤ t ≤ 2 s. To record the responded wave of the blocks with different seismic input methods, two measured points are set as shown in Figure 5a. Here, a₁, v₁, and d₁ represent the responsive horizontal acceleration, velocity, and displacement, respectively, of measured point 1; a₂, v₂, and d₂ represent the responsive horizontal acceleration, velocity, and displacement, respectively, of measured point 2.

As shown in Figure 5b, the multi-blocks Newmark method is used and the acceleration time history is applied to both block A and block B. It is appointed that right is positive and left is negative. The theoretical solution to this case is:

$$\begin{array}{cc}\hfill 1.\mathrm{If} & v_2=0\hfill \\ & \mathrm{If} a-g\tan \varphi >0, a_2=a-g\tan \varphi \hfill \\ & \mathrm{If} a+g\tan \varphi <0, a_2=a+g\tan \varphi \hfill \\ & \mathrm{If} a+g\tan \varphi \ge 0\ge a-g\tan \varphi , a_2=0\hfill \\ \hfill 2.\mathrm{If} & v_2>0, a_2=a-g\tan \varphi \hfill \\ \hfill 3.\mathrm{If} & v_2<0, a_2=a+g\tan \varphi \hfill \end{array}$$

(37)

where a is the input acceleration, g is the gravity acceleration, and φ is the friction angle of the contact face.

In Figure 5c, the large mass method is used, the horizontal constraint of block A is released, and acceleration of the input dynamic loading is applied to block A. In Figure 5d, the large stiffness method is used, and the input dynamic loading is applied to block A in the form of a horizontal displacement loading. In these cases, the responsive acceleration or displacement of block A is consistent with the input dynamic loading. The theoretical solution is:

$$\begin{array}{cc}\hfill 1.\mathrm{If} & v_1=v_2\hfill \\ & \mathrm{If} \left|a_1\right|\le g\tan \varphi , a_2=a_1\hfill \\ & \mathrm{If} a_1>g\tan \varphi , a_2=g\tan \varphi \hfill \\ & \mathrm{If} a_1<-g\tan \varphi ,a_2=-g\tan \varphi \hfill \\ \hfill 2.\mathrm{If} & v_1>v_2, a_2=g\tan \varphi \hfill \\ \hfill 3.\mathrm{If} & v_1<v_2, a_2=-g\tan \varphi \hfill \end{array}$$

(38)

Integrating on Equations (37) and (38), the analytical solutions of both velocity and displacement are obtained.

The responded waves of measured point 2 are illustrated in Figure 7. It shows that the DDA simulation results fit well with the corresponding analytical solutions, which indicates that these seismic input methods are correctly programmed in the seismic DDA method. In fact, for the multi-blocks Newmark method, it does not consider the seismic wave propagation and the seismic loads are improperly applied when the block falls. While for the large mass method and the large stiffness method, they consider the seismic wave propagation by timely passing the seismic loading input from the bottom block to the block system, and they are basically equivalent; since the bottom block is seen as a rigid foundation, these two methods are more suitable for simulating the seismic response of the geotechnical engineering with a rigid foundation, such as the shaking table test and the slope with clear bedrock [59].

5. Simulation of a Shaking Table Test Using the Seismic DDA Method

5.1. A Brief of the Shaking Table Test

The shaking table test is based on the large rock cavern complex of the Dagangshan Hydropower Station in Dadu River Basin, China [60]. According to the engineering layout of rock cavern complex, the physical model contains three main structures, including a main machine building, a main transformer chamber, and a tail surge chamber. The similarity ratio between the model and the prototype is 1:12.25. The structures and sizes of the physical model are shown in Figure 8. The similar materials of rock mass in the test are composed of iron powder, barite powder, quartz sand, gypsum, and water. The mass ratio of the five components is 176: 264: 66: 50: 60. To simulate the flexible boundary of the physical model, a polystyrene foam board is added inside the rigid model box [61].

During the shaking table test, a seismic wave measured in the Kobe earthquake of Japan which occurred in 1995 is used as the input load. According to the dynamic similarity relationship, the similarity ratio of the test time is 12.25. The upper limit of the effective frequency of the shaking table used is 50 Hz. Therefore, the waveform needs to be compressed by 12.25 times on the time axis, and the frequency components above 50 Hz are filtered out.

5.2. Numerical Simulation

According to the sizes of the physical model, the numerical model is established by using the DC module of the open-source program of 2D DDA. As shown in Figure 9, the DDA model includes rock mass, polystyrene foam board, rigid model box, shaking table, and a large mass block used for seismic input. There are 4381 blocks formed in the numerical model.

Table 3. Physical and mechanical parameters of various materials.

Materials	Density/(kg/m3)	Young’s Modulus/MPa	Poisson’s Ratio
Rock mass	2800	172	0.25
Polystyrene foam board	15	4.13	0.07
Steel	7800	21,000	0.30

lists the physical and mechanical parameters of various materials, where the rigid model box, shaking table, and large mass block are steel. The setting method of damping parameters is very important for DDA to simulate seismic response of rock mass engineering. The shaking table test results show that the damping ratio of similar materials of rock mass is about 0.04–0.05. In the DDA method, the damping ratio in the motion equation comes from numerical damping and viscous damping. The damping ratio of DDA is calculated as follows [62]:

$${\xi }_{total}={\xi }_{numerical}+{\xi }_{viscous}=\frac{\ln \left(1+\frac{{\Omega }^2}{2}\right)}{2\Omega }+\frac{\left(1-g{\mathrm{g}}^2\right)}{4\Omega }$$

(39)

where ξ is the damping ratio, $\Omega =\omega \Delta t$, gg is a dynamic coefficient, and ω is the circular frequency. In this paper, the time step is 0.004 s and the dynamic coefficient is set to 0.9965 for the rock mass. According to Equation (39), the distribution of damping ratio of the rock mass in frequency domain is illustrated. As shown in Figure 10, in the frequency domain of main concern (10–40 Hz), the DDA damping ratio is in good agreement with the shaking table test results. Similarly, the dynamic coefficients of polystyrene foam board and rigid model box are 0.9850 and 0.9995, respectively, and the damping ratios are controlled in the range of 10~20% and 1~2%, respectively.

During the DDA simulation, the large mass method is used to input the seismic wave. As shown in Figure 9, a large mass block at the bottom is modeled. The acceleration time history measured on the shaking table during the test is applied to the large mass block. At the same time, it must be ensured that the large mass block has no viscous damping, because damping will reduce the input seismic energy, and the zero damping for seismic input can be realized by setting the dynamic coefficient of the large mass block to 1.0.

5.3. Results Analysis

5.3.1. Seismic Input Method

In order to verify the correctness of the large mass method, as shown in Figure 9, a measured point at the large mass block (numbered 1) and a measured point at the shaking table (numbered 2) are set. The responded waves of measured points and the input seismic wave are illustrated in Figure 11. It shows that the responded waves are basically consistent with the input seismic wave. Therefore, the large mass method can effectively imitate the excitation principle of the shaking table, that is, the excitation system drives the vibration of the shaking table, and the shaking table drives the vibrations of the rigid model box and the rock mass.

5.3.2. Boundary Condition

In order to verify the effectiveness of polystyrene foam board as a flexible boundary, as shown in Figure 9, the responded acceleration of measured point 3 and measured point 4 in DDA model are recorded. As shown in Figure 12, the acceleration curves of the two measured points are basically the same. This means that the physical and mechanical parameters of the polystyrene foam board in numerical simulation are appropriate, which simulates the effect of the flexible boundary well.

5.3.3. Comparisons of Responded Acceleration of Rock Mass

In order to further verify the correctness of DDA simulation, As shown in Figure 9, taking the typical measured points of the main machine building structure as an example, including upstream side wall (measured point 5), downstream side wall (measured point 6), bottom plate (measured point 7), and crown arch (measured point 8), the responded accelerations of the test data and simulation results are compared. As shown in Figure 13, the test data and simulation results are basically consistent regardless of time history or Fourier spectrum. This means that the simulation results of DDA are reasonable.

5.3.4. Analysis of Propagation Law of Seismic Wave

In Figure 13, there is little difference between the responded accelerations of upstream side wall and downstream side wall, while the responded accelerations of bottom plate and crown arch are very different. This means that the propagation law of seismic wave has a strong correlation with the elevation. Based on this understanding, as shown in Figure 9, a measured line is arranged from the surface of the shaking table to the top of the model. During DDA simulation, the responded accelerations of 26 measured points with different elevations along the measured line are recorded. Here, it is assumed that the elevation of the shaking table surface is zero. In order to study the propagation law of the seismic wave, an amplification factor is defined as the ratio of the positive or negative maximum of responded acceleration to the positive or negative maximum of input seismic wave.

As shown in Figure 14, regardless of the positive or negative amplification factor, the amplification effect of seismic wave becomes more obvious with the increase of elevation, and the negative amplification factor is significantly greater than the positive amplification factor. According to the elevation, the distribution characteristics of magnification factor can be divided into two sections. The first section is the elevation less than 0.435 m, that is, below the bottom plate of the main machine building, and the amplification factor is in the range of 1.02–1.22, which indicates that the amplification effect is not obvious. The maximum positive amplification factor is 1.14, while the maximum negative amplification factor is 1.22. The second section is the elevation greater than 0.925 m, that is, above the crown arch of the main machine building, the positive amplification factor is in the range of 1.51–1.89, and the negative amplification factor is in the range of 2.31–3.20, which shows that the amplification effect is more obvious than that of the first section. It should be noted that during crossing the main machine building, due to the existence of a huge cavity, the seismic wave is reflected, refracted, and diffracted, and the amplification coefficient suddenly increases greatly.

6. Conclusions

The DDA method adopts the Newmark integration scheme to simulate the discontinuous motion process of block system in real time-steps. The small displacement assumption is satisfied in each time step, and the large displacement is formed by the accumulation of small displacement. This paper studies the solution principle of DDA method and the development status of a 2D open-source program. The seismic wave pre-processing and seismic input method are systematically investigated. The following conclusions are drawn:

(1) Based on the Newmark integration scheme, the integration algorithms of synthetic or measured seismic wave time history, correction function of seismic wave, and DDA simulation are unified.

(2) The formulae of various seismic input methods, including the multi-blocks Newmark method, large mass method, and large stiffness method, are deduced and implemented in DDA program using the C programming language. The large mass method and the large stiffness method are more suitable for rigid foundation, such as shaking table test.

(3) To simulate the shaking table test of a typical tunnel, the correctness of DDA in simulation of seismic input, boundary condition, and damping ratio of material are verified. At the same time, the propagation law of seismic wave in rock mass is analyzed. It reveals that the propagation of the seismic wave presents a significant amplification effect due to the reflection, refraction, and diffraction in the tunnel.