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

: 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 uniﬁed. 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 veriﬁed; 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 signiﬁcant ampliﬁcation effect due to the reﬂection, 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.


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 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. Expanding the above equation, the dynamic ordinary differential equations of the block system can be obtained, [M] .. 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: [M] ..
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: ..
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: where,

D n
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, where {∆D i } is the solution displacement of the block; (u 0 , v 0 ) is the block centroid for the rigid body translation; r 0 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: where ∆u i and ∆v i are components of the incremental displacements in the x direction and the y direction, respectively; Supposing there are n blocks in the defined block system, based on the Equation (7), the equilibrium equations of are written as:

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

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

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. 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, (a0, a1, a2,…, am−1, am). The time-step of the data is Δt, the total duration is T, the initial seismic velocity is V0 = 0, and the initial displacement d0 = 0. Using the Newmark integration scheme, the velocity can be obtained by integrating the acceleration once: The displacement can be obtained by integrating the acceleration twice: At the end of seismic action, the acceleration am may be non-zero, or vm and dm 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 It is assumed that the acceleration data observed by seismic stations or synthesized are known, (a 0 , a 1 , a 2 , . . . , 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 = 0, and the initial displacement d 0 = 0. Using the Newmark integration scheme, the velocity can be obtained by integrating the acceleration once: v n+1 = v n + ∆t (1 − γ)a n +γa n+1 (n = 0, 1, . . . , m − 1) The displacement can be obtained by integrating the acceleration twice: (1 − 2β)a n +2βa n+1 (n = 0, 1, . . . , m − 1) 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: 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: Equations (13)-(15) automatically satisfy the first condition (a). According the second condition (b): 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: 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: Using the relationship between the velocity and the acceleration described by Newmark integration scheme, Substituting Equation (19) into Equation (18), the velocity correction function is written as: Step 2: The displacement correction function is obtained from the correction function of acceleration and velocity: Using the relationship between the displacement and the acceleration described by Newmark integration scheme, Substituting Equation (22) into Equation (21), the displacement correction function is written as: 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.

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

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: 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: where 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:

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: [M] ..
where [M bb ] is the mass matrix of the bottom block, .. D g is the seismic acceleration vector. It can be decomposed into another form: 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): Both sides of Equation (29) are pre-multiplied by [M bb ] −1 and it becomes: ..
As a result of the large mass of the bottom block, diagonal elements of [M bb ] −1 tend to zero, if the damping of the bottom block is independent of the mass, Equation (29) can be approximately written as: ..
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.

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: [M] ..
where [K bb ] is the stiffness matrix of the bottom block, D g is the seismic displacement vector. It can be decomposed into another form: Expanding the second row of Equation (33): ..
Both sides of Equation (34) are pre-multiplied by [K bb ] −1 and it becomes: As a result of the rigid springs of the bottom block, diagonal elements of [K bb ] −1 tend to zero, if the damping of the bottom block is independent of the stiffness, Equation (35) can be approximately written as: 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.
As a result of the rigid springs of the bottom block, diagonal elements of [Kbb] −1 tend to zero, if the damping of the bottom block is independent of the stiffness, Equation (35) can be approximately written as: 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.

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 3 , and a gravity acceleration of 10 m/s 2 . 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.

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 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, a1, v1, and d1 represent the responsive horizontal acceleration, velocity, and displacement, respectively, of measured point 1; a2, v2, and d2 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: 1. If 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 1 , v 1 , and d 1 represent the responsive horizontal acceleration, velocity, and displacement, respectively, of measured point 1; a 2 , v 2 , and d 2 represent the responsive horizontal acceleration, velocity, and displacement, respectively, of measured point 2. As shown in Figure 6, the input dynamic loading is horizontal and its accelerat time history is a sine wave of a period a(t) = 10 sin(πt), where 0 ≤ t ≤ 2 s. To record responded wave of the blocks with different seismic input methods, two measu points are set as shown in Figure 5a. Here, a1, v1, and d1 represent the responsive h zontal acceleration, velocity, and displacement, respectively, of measured point 1; a2, and d2 represent the responsive horizontal acceleration, velocity, and displacement, spectively, of measured point 2. As shown in Figure 5b, the multi-blocks Newmark method is used and the acce ation time history is applied to both block A and block B. It is appointed that righ positive and left is negative. The theoretical solution to this case is: 1. If 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: 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: 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].
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: 1. If 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].  Table Test Using the Seismic DDA Method

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

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 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]: where ξ is the damping ratio, t  =   , gg is a dynamic coefficient, and ω is the circular 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.

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 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]: where ξ is the damping ratio, Ω = ω∆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       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 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.

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.

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. Figure 11. Input seismic wave and verification of large mass method.

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.

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.

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.

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.

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. ple, 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.  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. Figure 14. Amplification factor of responded acceleration at each measured point along the elevation.

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:

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.
Author Contributions: Conceptualization, X.F.; methodology, X.F. and Q.S.; validation, L.Z. and W.D.; formal analysis, J.K. and H.D.; writing-original draft preparation, X.F.; supervision, X.F. and Q.S.; funding acquisition, X.F. and Q.S. All authors have read and agreed to the published version of the manuscript.