*3.2. Implementation of Oblique Incident Waves*

#### 3.2.1. SV Waves

Based on the artificial boundary, the other key point in the numerical simulation of a dynamic problem is the implementation of seismic waves in finite element models. An efficient and high-precision method is to convert the seismic records, that is, acceleration records, together with the coefficients of the artificial boundary into equivalent nodal forces. Figure 3 presents a plane study area that was subjected to oblique incident SV waves. As shown in Figure 3, each truncated boundary was subjected to three waves that involved one incident wave and two reflected waves. When node A, one of the joints on the left truncated boundary, were considered, the incident seismic waves of node A reached the boundary directly (incident SV waves) or through reflection after first reaching the ground surface (reflected SV waves and reflected P waves). The displacement and velocity of the incident SV waves are denoted as *us*(*t*) and . *us*(*t*), respectively.

**Figure 3.** Sketch of the seismic waves on truncated boundary by obliquely incident SV waves.

Consequently, the free-field motion and the associated stress (in terms of *x* and *y* components) can be induced at the given point A (*x*0, *y*0) as follows:

$$\begin{cases} u\_{1x}(t) = u\_s(t - \Delta t\_1)\cos(\alpha) - A\_1 u\_s(t - \Delta t\_2)\cos(\alpha) + A\_2 u\_p(t - \Delta t\_3)\sin(\beta) \\\ u\_{1y}(t) = -u\_s(t - \Delta t\_1)\sin(\alpha) - A\_1 u\_s(t - \Delta t\_2)\sin(\alpha) + A\_2 u\_p(t - \Delta t\_3)\cos(\beta) \\\ \dot{u}\_{1x}(t) = \dot{u}\_s(t - \Delta t\_1)\cos(\alpha) - A\_1 \dot{u}\_s(t - \Delta t\_2)\cos(\alpha) + A\_2 \dot{u}\_p(t - \Delta t\_3)\sin(\beta) \\\ \dot{u}\_{1y}(t) = -\dot{u}\_s(t - \Delta t\_1)\sin(\alpha) - A\_1 \dot{u}\_s(t - \Delta t\_2)\sin(\alpha) + A\_2 \dot{u}\_p(t - \Delta t\_3)\cos(\beta) \end{cases} (8)$$

$$\begin{cases} \begin{array}{l} \sigma\_{lx} = S\_1 \left( \dot{u}\_s (t - \Delta t\_1) + S\_2 A\_1 \dot{u}\_s (t - \Delta t\_2) \right) + S\_3 A\_2 \dot{u}\_p (t - \Delta t\_3) \\\ \sigma\_{ly} = S\_4 \left( \dot{u}\_s (t - \Delta t\_1) + S\_5 A\_1 \dot{u}\_s (t - \Delta t\_2) \right) + S\_6 A\_2 \dot{u}\_p (t - \Delta t\_3) \end{array} \tag{9}$$

where *A*1 and *A*2 are the amplitude ratios (denoted as the reflected SV waves and reflected P waves to the incident SV waves), Δ*t* is the delay time of waves from the wavefront at *t* = 0 to the boundary node (node A), *α* and *β* are the angles between the vertical direction and reflected SV waves, and the reflected P waves, respectively. *A*1, *A*2 and *β* can be expressed as follows:

$$\begin{cases} A\_1 = \frac{c\_s^{-2} \sin 2a \sin 2\beta - c\_p^{-2} \cos^2 2a}{c\_s^{-2} \sin 2a \sin 2\beta + c\_p^{-2} \cos^2 2a} \\\ A\_2 = \frac{2c\_P c\_s \sin 2a \cos 2a}{c\_s^{-2} \sin 2a \sin 2\beta + c\_p^{-2} \cos^2 2a} \\\ \beta = \arcsin\left(\frac{c\_p \sin a}{c\_s}\right) \end{cases} \tag{10}$$

where *cs* and *cp* are the wave velocities of P waves and SV waves, respectively. The delay times Δ*t*1, Δ*t*2, Δ*t*3, and variables *S*1 to *S*6 are all boundary-dependent, which can be depicted as follows:

$$\begin{cases} \Delta t\_1 = y \cos a / c\_s\\ \Delta t\_2 = \left(2L\_y - y\right) \cos a / c\_s\\ \Delta t\_3 = \left(L\_y - y\right) / \left(c\_p \cos \beta\right) + \left(L\_y - \left(L\_y - y\right) \tan a \tan \beta\right) \cos a / c\_s\\ S\_1 = \frac{G \sin 2a}{c\_s}, S\_2 = -1, \ S\_3 = \frac{\lambda + 2G \sin^2 \beta}{\frac{c\_p}{c\_p}},\\ S\_4 = \frac{G \cos 2a}{c\_s}, S\_5 = 1, \ S\_6 = -\frac{G \sin 2\beta}{\frac{G}{c\_p}} \end{cases} \tag{11}$$

On the bottom boundary, they became:

$$\begin{cases} \Delta t\_1 = \mathbf{x} \sin a / c\_s\\ \Delta t\_2 = \left(2L\_y + \mathbf{x} \tan a\right) \cos a / c\_s\\ \Delta t\_3 = L\_y / \left(c\_p \cos \beta\right) + \left(L\_y \cos a + \mathbf{x} \sin a - L\_y \tan \beta \sin a\right) / c\_s\\ S\_1 = \frac{G \cos 2a}{c\_s}, \ S\_2 = 1, \ S\_3 = -\frac{G \sin 2\beta}{c\_p},\\ S\_4 = \frac{G \sin 2a}{c\_s}, \ S\_5 = 1, \ S\_6 = \frac{\lambda + 2G \cos^2 \beta}{c\_p} \end{cases} \tag{12}$$

On the right boundary, the stresses were the same but in opposite directions to those on the left boundary, and for the displacement, the delay times were Δ*t*1, Δ*t*2, and Δ*t*3, each of them increased by *Lx* sin *α*/*cs* for the additional travel distance.

#### 3.2.2. P Waves

Figure 4 presents a plane study area that was subjected to oblique incident P waves. Similar to the oblique incident SV waves, each truncated boundary was subjected to three waves that involved one incident wave and two reflected waves, as shown in Figure 4. When node B, one of the joints on the left truncated boundary, were considered, the incident seismic waves of node B reached the boundary directly (incident P waves) or through reflection after first reaching the ground surface (reflected P waves and reflected SV waves). The displacement and velocity of the incident P waves are denoted as *up*(*t*) and . *up*(*t*), respectively.

**Figure 4.** Sketch of the seismic waves on truncated boundary by obliquely incident P waves.

Consequently, the free-field motion and associated stress (in terms of their *x* and *y* components) can be induced at the given point B (*x*1, *y*1) as follows:

$$\begin{cases} u\_{lx}(t) = u\_p(t - \Delta t\_4)\sin(a) + A\_3 u\_p(t - \Delta t\_5)\sin(a) + A\_4 u\_s(t - \Delta t\_6)\cos(\beta) \\\ u\_{ly}(t) = u\_p(t - \Delta t\_4)\cos(a) - A\_3 u\_p(t - \Delta t\_5)\cos(a) + A\_4 u\_s(t - \Delta t\_6)\sin(\beta) \\\ \dot{u}\_{lx}(t) = \dot{u}\_p(t - \Delta t\_4)\sin(a) + A\_3 \dot{u}\_p(t - \Delta t\_5)\sin(a) + A\_4 \dot{u}\_s(t - \Delta t\_6)\cos(\beta) \\\ \dot{u}\_{ly}(t) = \dot{u}\_p(t - \Delta t\_4)\cos(a) - A\_3 \dot{u}\_p(t - \Delta t\_5)\cos(a) + A\_4 \dot{u}\_s(t - \Delta t\_6)\sin(\beta) \end{cases} (13)$$

$$\begin{cases} \begin{array}{l} \sigma\_{lx} = S\_7 \Big( \dot{u}\_p(t - \Delta t\_4) + S\_8 A\_3 \dot{u}\_p(t - \Delta t\_5) \Big) + S\_9 A\_4 \dot{u}\_s(t - \Delta t\_6) \\\ \sigma\_{ly} = S\_{10} \Big( \dot{u}\_p(t - \Delta t\_4) + S\_{11} A\_3 \dot{u}\_p(t - \Delta t\_5) \Big) + S\_{12} A\_4 \dot{u}\_s(t - \Delta t\_6) \end{array} \end{cases} \tag{14}$$

where *A*3 and *A*4 are the amplitude ratios (denoted as the reflected P waves and reflected SV waves to the incident P waves). Δ*t* denotes the delay time of waves from the wavefront at *t* = 0 to the boundary node (node B). *α* and *β* are the angles between the vertical direction and reflected P waves together with the reflected SV waves, respectively. *A*3, *A*4 and *β* can be expressed as follows:

$$\begin{cases} A\_3 = \frac{c\_s^{-2} \sin 2a \sin 2\beta - c\_p^{-2} \cos^2 2\beta}{c\_s^{-2} \sin 2a \sin 2\beta + c\_p^{-2} \cos^2 2\beta} \\\ A\_4 = \frac{2c\_p c\_s \sin 2a \cos 2\beta}{c\_s^{-2} \sin 2a \sin 2\beta + c\_p^{-2} \cos^2 2\beta} \\\ \beta = \arcsin\left(\frac{c\_s \sin a}{c\_p}\right) \end{cases} \tag{15}$$

where *cs* and *cp* are the wave velocities of P waves and SV waves, respectively. The delay times Δ*t*4, Δ*t*5, Δ*t*6, and the variables *S*7 to *S*12 are all boundary-dependent, which can be depicted as follows:

$$\begin{cases} \Delta t\_4 = y \cos a / c\_p \\ \Delta t\_5 = \left( 2L\_y - y \right) \cos a / c\_p \\ \Delta t\_6 = \left( L\_y - y \right) / \left( c\_s \cos \beta \right) + \left( L\_y - \left( L\_y - y \right) \tan a \tan \beta \right) \cos a / c\_p \\ \end{cases} \tag{16}$$

$$\begin{array}{l} S\_7 = \frac{\lambda + 2G \sin^2 a}{c\_p}, \ S\_8 = 1, \ S\_9 = \frac{G \sin 2\beta}{c\_s}, \\\ S\_{10} = \frac{G \sin 2a}{c\_p}, \ S\_{11} = -1, \ S\_{12} = -\frac{G \cos 2\beta}{c\_s} \end{array} \tag{17}$$

On the bottom boundary, they became:

$$\begin{cases} \Delta t\_4 = x \sin a / c\_p \\ \Delta t\_5 = \left( 2L\_y + x \tan a \right) \cos a / c\_p \\ \Delta t\_6 = L\_y / (c\_s \cos \beta) + \left( L\_y \cos a + x \sin a - L\_y \tan \beta \sin a \right) / c\_p \\ S\tau = \frac{G \sin 2a}{c\_p}, \ S s = -1, \ S s = -\frac{G \cos 2\beta}{c\_s}, \\ S\_{10} = \frac{\lambda + 2G \cos^2 a}{c\_p}, \ S\_{11} = 1, \ S\_{12} = -\frac{G \sin 2\beta}{c\_s} \end{cases} \tag{17}$$

On the right boundary, the stresses were the same but in directions opposite to those on the left boundary, and for the displacement, the delay times were Δ*t*4, Δ*t*5, and Δ*t*6" each of which increased by *Lx* sin *α*/*cp* for the additional travel distance.

#### *3.3. Verification*

The incident seismic waves were implemented using the ABAQUS software, together with a self-developed MATLAB program. The specific implementation is described in the flowchart illustrated in Figure 5. Two numerical test examples were used to verify the validity of the oblique incident seismic waves. One was the propagating process of oblique incident waves in a truncated region, which was considered to explore the accuracy of the input method in a semi-infinite field. The other was the response of a semicircular canyon input by oblique incident waves, which was aimed at determining the accuracy of the simulation in the topography.

**Figure 5.** Flow chart for input seismic waves implemented in ABAQUS.

#### 3.3.1. Test Example 1

A truncated rectangular domain was adopted to simulate the propagation process of oblique incident waves in a semi-infinite field, as depicted in Figure 6a. The sizes of the truncated regions were as follows: *Lx* = 800 m and *Ly* = 400 m. The region was sufficiently large to detect the propagation of all the incident and reflected waves. The incident angles of the P waves and SV waves were assumed to be 30◦ and 18◦, respectively.

**Figure 6.** Numerical model and incident waves of test example 1: (**a**) finite element model and (**b**) displacement time history of incident waves.

The material parameters of the domain were as follows: Poisson's ratio = 0.25, elastic modulus = 1.25 Gpa and mass density, *ρ* = 2000 kg/m3. Points A (400, 400) and B (400, 0) are the observation points. An impulse was adopted as an incident seismic wave, as shown in Figure 6b. The impulse equation is defined as follows [57]:

$$P(t) = 16P\_0\left[G(t) - 4G\left(t - \frac{1}{4}\right) + 6G\left(t - \frac{1}{2}\right) - 4G\left(t - \frac{3}{4}\right) + G(t - 1)\right] \tag{18}$$

where *G* = (*t*/*T*0) 3*H*(*t*/*T*0) *H*(*t*) is the Heaviside function, *P*0 is the peak value of an impulse, *P*0 is set to 1.0 m; and *T*0 is the acting time of the impulse, where *T*0 = 0.25 s.

Figure 7 shows the displacement contours with the time of incident SV waves. As illustrated in the figures, the propagation processes of the incident waves and reflected waves at the ground surface were simulated effectively in a semi-infinite field. In addition,

the displacement time histories of the observation points (points A and B) are shown in Figure 8. As shown in the figures, the numerical results proposed by the input method were in good agreement with the analytical solutions, where the analytical solutions were calculated using the elastic wave propagation theory [58,59].

**Figure 8.** Displacement time history at observation points of incident SV waves: (**a**) horizontal displacement, and (**b**) vertical displacement (solid lines and dashed lines are numerical results of points A and B, respectively, with circles showing analytical solutions).

The displacement contours with the time of P waves in Figure 9 show the propagation processes of the incident waves and reflected waves at the ground surface. Meanwhile, the displacement contours and the displacement time histories of the observation points (points A and B) in Figure 10 both confirm that the input method was effective and precise in the semi-infinite field. Therefore, the oblique incident SV waves and P waves in a semi-infinite field could be simulated using the proposed input method, and the ability of the viscous-spring artificial boundary to absorb the scattered waves was also verified.

**Figure 9.** Displacement contours of P waves propagated in semi-infinite space.

**Figure 10.** Displacement time history at observation points of incident P waves: (**a**) horizontal displacement, and (**b**) vertical displacement (solid lines and dashed lines are numerical results of points A and B, respectively, with circles showing analytical solutions).

3.3.2. Test Example 2

A truncated rectangular domain with a semicircular canyon is shown in Figure 11a. The truncated region was 0 ≤ *x* ≤ 200 m, 0 ≤ *y* ≤ 200 m and the radius of the canyon was 25 m. The material parameters of the ground were as follows: Poisson's ratio = 0.3, elastic modulus = 208 MPa and mass density, *ρ* = 2000 kg/m3. A Ricker wavelet was adopted as the incident wave. The acceleration-time history of the wave is illustrated in Figure 11b. The time history of the Ricker wavelet was defined as follows:

$$R(t) = \left(1 - 2\pi^2 f\_0^2 t^2\right) \exp\left(-\pi^2 f\_0^2 t^2\right) \tag{19}$$

with an excitation frequency *f* 0 of 4.0 Hz.

The incident angles had values of 25◦ under SV waves and 30◦ under P waves. Figures 12 and 13 show the horizontal displacement amplification ratios Uh and vertical displacement amplification ratios Uv of the ground surface on the canyon with incident SV waves and incident P waves, respectively. Uh is defined as the ratio of the horizontal displacement on the ground surface to the input horizontal displacement, and Uv is the ratio of the vertical displacement. The numerical results presented in Figures 12 and 13 were in good agreement with the analytical solutions proposed by Wong [60]. Therefore, in the topographic analysis, the proposed input method could effectively simulate the propagation process of oblique incident seismic waves, and the viscous-spring artificial boundary could absorb the scattered waves perfectly because of the topography.

**Figure 11.** Numerical model and input waves of test example 2: (**a**) finite element model and (**b**) acceleration time-history of incident waves.

**Figure 12.** *Cont*.

**Figure 12.** Displacement amplification ratios of ground under incident SV-waves with angle of 25◦: (**a**) horizontal ratios and (**b**) vertical ratios (solid lines are analytical solutions by Wong [60] and circles are numerical results of finite element method based on ABAQUS).

**Figure 13.** Displacement amplification ratios of ground under incident P-waves with angle of 30◦: (**a**) horizontal ratios and (**b**) vertical ratios (solid lines are analytical solutions by Wong [60] and circles are numerical results of finite element method based on ABAQUS).