## 2 Numerical Modeling Using 3D Discrete Elements
### Setting the Dynamic Parameters  
Unlike the general static case, when numerical analysis is performed under dynamic conditions careful consideration must be given to the relevant dynamic parameters and conditions so that appropriate values can be chosen. Usually, the model must first be run under static (i.e., gravity-only) conditions until static equilibrium is reached. Thereafter, the dynamic analysis can be performed. The grid size of the model is set to 2 m in order to satisfy the relationship \(\Delta l\leq(1/10\sim 1/8)\gamma\), which is necessary to ensure the accurate propagation of the stress waves within the rock mass [21], where \(\Delta l\)--grid size and \(\lambda\)--wavelength. Rayleigh damping is also introduced into the dynamic analysis to account for damping phenomena:  
\[[C]=\xi_{\text{min}}\omega_{\text{min}}[M]+\xi_{\text{min}}[K]/\omega_{\text{ min}} \tag{3}\]  
where \([M]\) is the mass matrix, \([K]\) is the stiffness matrix, \(\xi_{\text{min}}\) is the critical damping ratio, and \(\omega_{\text{min}}\) is the circular frequency. In this paper, the values \(\xi_{\text{min}}=0.02\) and \(\omega_{\text{min}}=5\) Hz are adopted. The bottom of the model was chosen to be a viscous (no-reflection) boundary, that is, separate dampers were set in the normal and tangential directions at the bottom of the model to provide viscous normal and tangential traction, respectively. Free field boundaries were adopted on the four sides of the model to simulate free field motion (Figure 3). \begin{table}
\begin{tabular}{c c c c c c c c c} \hline \hline  & **Density** & **Young’s** & **Poisson’s** & **Cohesion** & **Friction** & **Tensile** & **Dilation** & **Normal** & **Tangential** \\  & **(kN/m3)** & **(GPa)** & **Ratio** & **(MPa)** & **Angle (°)** & **Strength** & **Angle (°)** & **Angle (°)** & **St stiffness** \\ \hline Rock mass & 25.97 & 10.5 & 0.25 & 1.5 & 45 & 0.5 & 0 & – & \\ Structural planes & – & – & – & 0.24 & 32 & 0.001 & 0 & 58.33 & 19.44 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Physical and mechanical parameters used for the rock mass and structural planes. Figure 3: Schematic diagrams illustrating the dynamic boundary conditions imposed on the 3D discrete-element models showing: (**a**) a side view of the model, and (**b**) a 3D spatial view. As the bottom of the model is set to be a viscous boundary, the dynamic load should be converted from acceleration history to velocity history (by integration) and then to stress history according to:  
\[\sigma_{n} =2(\rho C_{p})v_{n} \tag{4}\] \[\tau_{s} =2(\rho C_{s})v_{s} \tag{5}\]  
which are applied to the bottom of the model, where \(\sigma_{n}\) and \(\tau_{s}\) are the applied normal and shear stresses, respectively, \(\rho\) is the rock density, \(C_{p}\) and \(C_{s}\) are the speeds of the P- and S-waves, respectively, and \(v_{n}\) and \(v_{s}\) are the normal and tangential velocities, respectively. In this paper, two types of dynamic loads were used in the dynamic analysis corresponding to sinusoidal and natural seismic waves. The whole dynamic process lasted 6 s. The expressions for the acceleration history of the sinusoidal waves in the three directions are all in the form: \(a_{t}=A\sin(2\pi ft)\)--see Figure 4--the corresponding velocity history is therefore given by \(v_{t}=A\cos(2\pi ft)/2\pi f\)). In this paper, the parameter values adopted are:\(f=5\) Hz for the frequency, because the main energy of the earthquake is concentrated in this frequency band; and \(A=3.14\) m/s\({}^{2}\) for the amplitude as these values correspond to a natural Wechuan earthquake in 2008, of intensity \(M_{s}=8.0\)[22].