Soil seismic analysis for 2D oblique incident waves using exact free-field responses by frequency-based finite/infinite element method

Abstract

The seismic analysis of a viscoelastic half-space under two-dimensional (2D) oblique incident waves is carried out by the finite/infinite element method (FIEM). First, the frequency-domain exact solutions for the displacements and stresses of the free field are derived in general form for arbitrary incident P and SV waves. With the present formulation, no distinction needs to be made for SV waves with over-critical incident angles that make the reflected P waves disappear, while no critical angle exists for P waves. Next, the equivalent seismic forces of the earthquake (Taft Earthquake 1952) imposed on the near-field boundary are generated by combining the solutions for unit ground accelerations with the earthquake spectrum. Based on the asymmetric finite/infinite element model, the frequency-domain motion equations for seismic analysis are presented with the key parameters selected. The results obtained in frequency and time domain are verified against those of Wolf’s, Luco and de Barros’ and for inversely computed ground motions. The parametric study indicated that distinct phase difference exists between the horizontal and vertical responses for SV waves with over-critical incident angles, but not for under-critical incident angles. Other observations were also made for the numerical results inside the text.

Graphical abstract

Keywords

oblique incident waves / critical angle / half-space / finite/infinite element approach / seismic response analysis

1 Introduction

The dynamic response of a site hit by seismic waves with oblique incident angles is significantly different from those by vertical incident seismic waves. The irregular response owing to oblique incident waves cannot be ignored, which may seriously threaten the normal operation and life of the affected area. It was pointed out that oblique incident seismic waves often occur in the half-space with hard soils, whose effect on structural responses should be considered in the seismic design [1]. Previously, various theoretical and numerical approaches have been presented to study the half-space response to the seismic waves with arbitrary incident angles.

As for theoretical aspect, the elastic wave behavior in the half-space with various forms of layer medium was studied by Ewing et al. [2], and in elastic solids by Graff [3]. Wolf [4] solved the soil-structure interaction under incident seismic body waves systematically. Lee and Karl [5] used the Fourier-Bessel series to analyze the scattering and diffraction of plane SV waves in an elastic half-space with underground cavities embedded at different depths. Yuan and Liao [6] derived the closed-form solution of two-dimensional (2D) scattering of plane SH waves for a half-space containing the cylindrical alluvial valley by the wave functions expansion. Davis et al. [7] combined the Fourier-Bessel series with a convex approximation for the free surface of the half-space, to present the theoretical transverse response of unlined cavities in an elastic half-space. Using the same method, Liang et al. [8] derived the theoretical solution for the incidence of P waves scattering caused by a circular-arc canyon with a covering layer.

Due to the lack of versatility of the analytical methods, various numerical methods were called on to explore the behavior of wave propagation in half-space with variations in local topography or underground structures. In this concern, the half-space is divided into two parts as the near field, containing the underground structures and cavities of interest, and the far field that extends to infinity. Previously, a generalized inverse method for approximating the boundary conditions was proposed by Wong [9] to consider the diffraction of P, SV and Rayleigh waves and wave propagation near the semi-elliptical canyon. Wong et al. [10] analyzed the dynamic stress and displacement response of a cylindrical tunnel embedded in an elastic half-space by using the eigenfunction expansions coupled with the finite element method (FEM). de Barros and Luco applied an indirect boundary integral method to examining the diffraction of oblique incident waves in a half-space induced by a cylindrical cavity [11] and the response for an infinitely long cylindrical cavity [12]. Then they performed the seismic analysis of a layered half-space with an embedded cylindrical shell [13,14] and expanded it into the three-dimensional (3D) form [15]. Subsequently, Stamos and Beskos used the boundary element method (BEM) to analyze the dynamic response of large 3D underground structure [16] and the 3D seismic response of lined tunnels with uniform characteristics along the direction perpendicular to the cross section in the half-space [17]. Until now, the BEM-based methods have been widely used by scholars in dealing with various half-space problems, including the analysis of shear deformable plates on elastic half-space, layered or saturated soils, valley or hill, 2D or 3D waves, see Shaaban and Rashed [18], Ba et al. [19–22] and Liang et al. [23], as well as the seismic ground response of the half-space with the subsurface box-shaped lined tunnel [24]. Recently, Mostafa Shaaban et al. [25–29] studied the wave propagation in the frequency domain for acoustic problems coupled with the boundary integral equation. Although the BEM has high accuracy in computation, the selection of proper Green’s functions depends on the specific problem to be considered.

In the literature, various methods and considerations for soil dynamic analysis are available. Komatitsch et al. [30] presented the spectral element method (SEM) to simulate the elastic wave propagation behavior in 2D and 3D realistic geological structures. Liu et al. [31,32] developed the time domain based viscous-spring artificial boundary for treating the soil-structure interaction. Based on the scaled boundary finite element method (SBFEM), Du and Lin [33] proposed an improved time domain numerical method to analyze the structure-foundation interaction problem. Focusing on rational approximation of the frequency response function, Du and Zhao [34] proposed an effective method to analyze the temporally local representation of unbounded soil. Hatzigeorgiou and Beskos [35] conducted the 3D seismic inelastic analysis of lined tunnels considering the soil-structure interaction by combining the viscous absorbing boundary with the FEM. Yu et al. [36] presented a multiscale method to analyze long tunnels by combining the coarse- and fine-scale finite elements, which allows more detailed response in locations of potential damage or of interest to be obtained. Zhao et al. [37] proposed a one-dimensional (1D) finite element artificial boundary method for simulating the homogeneous plane elastic wave propagation in a layered half-space caused by the plane waves with oblique incident angles. Besides, scholars also used the exponential basis functions to construct and apply the time-domain based absorbing boundary conditions to the wave propagation problems [38,39].

In previous soil dynamic analyses, viscous-spring artificial boundary has often been adopted. For instance, it was used by Huang et al. [40] in studying the impact of incident angles of earthquakes on the seismic response of long lined tunnels, and by Yan et al. [41] in studying the transverse seismic response of lined circular tunnels caused by asynchronous oblique P and SV waves. This is a convenient approach, but the physical properties of springs and dampers to be imposed on the artificial boundary require further justification.

Infinite elements have also been used for simulating the infinite boundary, for their relative ease in implementation, and for the fact that the parameters involved can be given in a more rational way. Ungless [42] and Bettess [43] proposed the concept of infinite elements for solving problems with infinity spatial domain. Zhao et al. [44] used the finite/infinite element method (FIEM) to simulate the scattering of incident P- and SV-waves in infinite media, and those of canyon topographic and geologic conditions [45,46], along with an effective wave input procedure for ground motions [47]. Yang et al. [48] combined the finite with the infinite elements to analyze the soil vibration by line loads, in which the key parameters for analysis using the finite and infinite elements were clearly given, such as the required element size, mesh range, and selection of wavenumbers and attenuation coefficients. Yun and co-workers derived various frequency-domain infinite elements to simulate different parts of the far field of a 2D layered half-space with soil medium [49], and then extended this approach to the time domain analysis [50].

Concerning the FIEM, Wang et al. [51] used the FIEM to study the vertical vibrations in unbounded saturated composite foundation with the axisymmetrical infinite elements proposed. Kouroussis et al. [52] studied the generation and propagation of ground vibrations under the action of railway traffic based on the vehicle-track-soil model. Yang et al. [53] used it to study the 2D response of an elastic half-space containing a circular cavity under vertical harmonic P and SV waves, focused on comparison of the results with previous solutions. For its versatility, the method was also employed to analyze the seismic response of single and double tunnels to vertical P and SV waves, respectively [54,55]. In the aforementioned studies using the FIEM, the seismic waves were all assumed to vertically hit the ground. As was mentioned previously, there are situations where the effect of oblique incident seismic waves cannot be ignored. Therefore, it is necessary to expand the FIEM such that oblique incident seismic waves can be accommodated, especially for SV waves with arbitrary incident angles. This requires a revision of the fundamental theory underlying the method of analysis.

To fill the above gap, the fundamental theory of a viscoelastic half-space subjected to 2D oblique incidence of seismic waves will be employed. The paper is organized as follows. In Section 2, the stresses and displacements of the free field are related to the Helmholtz potential for P- and SV-wave incidence, and for solution, the partial differential wave equations are transformed to the ordinary ones by the Fourier transformation. In Section 3, the exact solutions of the free-field displacements and stresses in frequency-space domain under either the oblique incident P or SV waves are derived. Particularly, the critical angle (above which the reflected P waves will disappear) for the incident SV waves is included with no extra effort. In Section 4, for a unit horizontal and vertical accelerations on the free surface (for P and SV waves, respectively), the displacements and stresses are computed for the near-field boundary using the above exact solutions. Then, the equivalent seismic forces acting on near-field boundary are calculated using Zhao and Valliappan’s formula [47] for the spectral frequencies considered. In Section 5, the finite/infinite element model under asymmetric seismic waves is established, the key parameters are selected, and the seismic analysis procedure in frequency domain is outlined. In Section 6, the validity of the proposed procedure is verified against Wolf’s and Luco and de Barros’ solutions, as well as for inversely computed ground motions. Finally, in Section 7, a parametric study for the seismic response of the half-space in frequency and time domains is carried out, with focus placed on the response for SV waves with over-critical incident angles. The last Section 8 is the conclusions.

2 Governing equations of motion for 2D half-space

A 2D (x-y) homogeneous and isotropic half-space, as shown in Fig.1, is considered for the seismic analysis under oblique incident P or SV waves. The equation of motion for an elastic solid in terms of the displacement field u is [56]:

(1)(\lambda +\mu )\mathrm{\nabla }\mathrm{\nabla }\cdot \mathit{u}+\mu {\mathrm{\nabla }}^2 \mathit{u}+\rho \mathit{f}=\rho {\mathit{u}}^{″},

where \lambda and \mu denote Lamé’s constants, and \rho the mass density and \mathit{f} the body force of the solid, respectively. The displacement of a generic point \mathit{u}=(u,v) consists of two components along the horizontal and vertical directions (x- and y-axes). Through use of the Helmholtz potential, Eq. (1) can be replaced by a set of equations that can be easily managed [56]. Namely, the propagation of the compression waves can be expressed by a scalar potential \Phi \left(\mathit{x},t\right), and that of the shear waves by a vector potential, \mathit{\Psi }\left(\mathit{x},t\right)=\left({\Psi }_{xy}^{\mathrm{S}\mathrm{V}},{\Psi }_{xy}^{\mathrm{S}\mathrm{H}}\right), with \mathit{x} denoting the coordinate vector of the generic point. The displacement field \mathit{u} can be represented as

\begin{array}{}\text{(2a)}& \mathit{u}(u,v)=\mathrm{\nabla }\Phi \left(\mathit{x},t\right)+\mathrm{\nabla }\times \mathit{\Psi }\left(\mathit{x},t\right),\end{array}

\begin{array}{}\text{(2b)}& {\Psi }_{xy}^{\text{SH}}=0.\end{array}

The condition {\Psi }_{xy}^{\mathrm{S}\mathrm{H}}=0 indicates that no out-of-plane motion is caused by the incident SH waves. Consequently, the two components of \mathit{u} can be uniquely determined from the potentials \Phi and {\Psi }_{xy}^{\mathrm{S}\mathrm{V}} for the 2D shear waves on the x-y plane. By neglecting the body force, one can substitute Eq. (2) into Eq. (1) to obtain the wave equations in partial differential form [56], i.e.,

\begin{array}{}\text{(3a)}& {\mathrm{\nabla }}^2\Phi \left(\mathit{x},t\right)=\frac{1}{c_{\text{P}}^2}\frac{{\mathrm{\partial }}^2\Phi \left(\mathit{x},t\right)}{\mathrm{\partial }{t}^2},\end{array}

\begin{array}{}\text{(3b)}& {\mathrm{\nabla }}^2\mathit{\Psi }\left(\mathit{x},t\right)=\frac{1}{c_{\text{S}}^2}\frac{{\mathrm{\partial }}^2\mathit{\Psi }\left(\mathit{x},t\right)}{\mathrm{\partial }{t}^2},\end{array}

in which the velocities of the P and SV waves propagating in the solid, c_{\mathrm{P}} and c_{\mathrm{S}}, respectively, are [56]

\begin{array}{}\text{(4a)}& c_{\text{P}}=\sqrt{\left(\lambda +2\mu \right)/\phantom{\left(\lambda +2\mu \right)\rho }\rho },\end{array}

\begin{array}{}\text{(4b)}& c_{\text{S}}=\sqrt{\mu /\phantom{\mu \rho }\rho }.\end{array}

In this study, the soil is regarded as a viscoelastic medium. According to Seed and Idriss [57], the Lamé’s constants, \lambda and \mu, are replaced by

\begin{array}{}\text{(5a)}& \lambda ={\lambda }_{\text{l}}\left[1+2\mathrm{i}\beta \mathrm{s}\mathrm{g}\mathrm{n}\left(\omega \right)\right],\end{array}

\begin{array}{}\text{(5b)}& \mu ={\mu }_{\text{l}}\left[1+2\mathrm{i}\beta \mathrm{s}\mathrm{g}\mathrm{n}\left(\omega \right)\right],\end{array}

in which {\lambda }_{\mathrm{l}} and {\mu }_{\mathrm{l}} denote the Lamé’s constants of the linearly elastic solid, i is the imaginary unit, \omega is the circular frequency (rad/s), and \beta is the hysteretic damping ratio.

The two displacement components (u,v) of the half-space can be obtained from Eq. (2) as:

\begin{array}{}\text{(6a)}& u=\frac{\mathrm{\partial }\Phi }{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\Psi }_{xy}^{\text{SV}}}{\mathrm{\partial }y},\end{array}

\begin{array}{}\text{(6b)}& v=\frac{\mathrm{\partial }\Phi }{\mathrm{\partial }y}-\frac{\mathrm{\partial }{\Psi }_{xy}^{\text{SV}}}{\mathrm{\partial }x}.\end{array}

Then, one can combine the strain-displacement relations with Hooke’s law to obtain the stresses of a viscoelastic half-space, in terms of Φ and {\Psi }_{xy}^{\mathrm{S}\mathrm{V}}, as:

\begin{array}{}\text{(7a)}& {\sigma }_{xx}=\lambda \left(\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}\right)\Phi +2\mu \left(\frac{{\mathrm{\partial }}^2\Phi }{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{\Psi }_{xy}^{\text{SV}}}{\mathrm{\partial }x\mathrm{\partial }y}\right),\end{array}

\begin{array}{}\text{(7b)}& {\sigma }_{yy}=\lambda \left(\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}\right)\Phi +2\mu \left(\frac{{\mathrm{\partial }}^2\Phi }{\mathrm{\partial }{y}^2}-\frac{{\mathrm{\partial }}^2{\Psi }_{xy}^{\text{SV}}}{\mathrm{\partial }x\mathrm{\partial }y}\right),\end{array}

\begin{array}{}\text{(7c)}& {\tau }_{xy}=2\mu \frac{{\mathrm{\partial }}^2\Phi }{\mathrm{\partial }x\mathrm{\partial }y}+\mu \left(\frac{{\mathrm{\partial }}^2{\Psi }_{xy}^{\text{SV}}}{\mathrm{\partial }{y}^2}-\frac{{\mathrm{\partial }}^2{\Psi }_{xy}^{\text{SV}}}{\mathrm{\partial }{x}^2}\right).\end{array}

The double Fourier transform {\hat{H}}_{\mathrm{n}} (denoted with “\hat{}”) and its inverse to be adopted for the wave function H_{\mathrm{n}} are defined as follows:

\begin{array}{rl}{\hat{H}}_{\text{n}}(k_x,y,\omega )=& \frac{1}{{(2\text{π})}^2}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{H}_{\text{n}}(x,y,t)\\ & \exp (\text{i}{k}_{x}x)\exp (-\text{i}\omega t)\text{d}x\text{d}t,\end{array}

\begin{array}{}\text{(8b)}& H_{\text{n}}(x,y,t)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\hat{H}}_{\text{n}}(k_x,y,\omega )\exp (-\text{i}{k}_{x}x)\exp (\text{i}\omega t)\text{d}{k}_x\text{d}\omega ,\end{array}

where k_x denotes the wave number along the x-axis. It is noted that the term \mathrm{e}\mathrm{x}\mathrm{p} (-\mathrm{i}{k}_{x}x ) in Eq. (8b) contains the minus sign in the exponent, but inversely in Eq. (8a). This condition ensures that the seismic waves will attenuate along the positive x-axis. With the double Fourier transformation performed, the wave equations in Eq. (3) can be converted into the following ordinary differential equations with the vertical coordinate y as one of the variables:

\begin{array}{}\text{(9a)}& \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}\hat{\Phi }(k_x,y,\omega )+\left(-k_x^2+k_{\mathrm{P}}^2\right)\hat{\Phi }(k_x,y,\omega )=0,\end{array}

\begin{array}{}\text{(9b)}& \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}{\hat{\Psi }}_{xy}^{\mathrm{S}\mathrm{V}}(k_x,y,\omega )+\left(-k_x^2+k_{\text{S}}^2\right){\hat{\Psi }}_{xy}^{\text{SV}}(k_x,y,\omega )=0,\end{array}

where k_{\mathrm{P}}=\omega /c_{\mathrm{P}}, k_{\mathrm{S}}=\omega /c_{\mathrm{S}}.

3 Free field under oblique incident seismic waves

3.1 Exact solution for incident P waves

Based on Snell’s law, both the phenomena of wave reflection and waveform conversion will occur on the ground under the action of P-wave incidence (P_{\mathrm{I}}) [2,3]. Namely, both the reflected P and SV waves (P_{\mathrm{R}} and {SV}_{\mathrm{R}}) are observed as shown in Fig.2, in which the incident and reflected angles of the P waves are denoted by {\theta }_{\mathrm{P}\mathrm{i}} and {\theta }_{\mathrm{P}\mathrm{P}\mathrm{r}}, respectively, and the reflected angle of the SV waves by {\theta }_{\mathrm{P}\mathrm{S}\mathrm{r}}. According to Snell’s law, the incident waves have the same wave number along the x-axis as that of the reflected ones for incident P waves [2]. For brevity, let us use k_{x\mathrm{P}} to denote the wave numbers of all body waves along the x-axis under oblique incident P waves, i.e., with k_{x\mathrm{P}}=\omega \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{P}\mathrm{i}} /c_{\mathrm{P}}.

By letting r_{\mathrm{P}}^2= k_{x\mathrm{P}}^2- k_{\mathrm{P}}^2, r_{\mathrm{S}}^2= k_{x\mathrm{P}}^2- k_{\mathrm{S}}^2, the solutions of {\hat{\Phi }}_{\mathrm{P}} and {\hat{\Psi }}_{xy,\mathrm{P}}^{\mathrm{S}\mathrm{V}} obtained from Eq. (9) are

\begin{array}{}\text{(10a)}& {\hat{\Phi }}_{\text{P}}=I_{\text{P}}\exp \left(r_{\text{P}}y\right)+R_{\text{P}}\exp \left(-r_{\text{P}}y\right),\end{array}

\begin{array}{}\text{(10b)}& {\hat{\Psi }}_{xy,\text{P}}^{\text{SV}}=R_{\text{S}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-r_{\text{S}}y\right),\end{array}

where I_{\mathrm{P}} represents the amplitude of the incident P waves and R_{\mathrm{P}} that of the reflected ones, and R_{\mathrm{S}} is the amplitude of the reflected SV waves. Meanwhile, the term with I_{\mathrm{P}} denotes the waves that propagate and attenuate upward, which are accounted for by \mathrm{e}\mathrm{x}\mathrm{p} \left(r_{\mathrm{P}}y\right), the term with R_{\mathrm{P}} the waves downward by \mathrm{e} \mathrm{x}\mathrm{p}\left({-r}_{\mathrm{P}}y\right), and the term with R_{\mathrm{S}} by \mathrm{e}\mathrm{x}\mathrm{p} \left({-r}_{\mathrm{S}}y\right).

Subsequently, the transformed displacements and stresses can be obtained by substituting the potential functions in Eq. (10) into Eqs. (6) and (7), and written in terms of I_{\mathrm{P}}, R_{\mathrm{P}}, and R_{\mathrm{S}} as follows:

\begin{array}{}\text{(11a)}& {\left\{\hat{u},\hat{v}\right\}}_{\text{P}}^{\text{T}}=\left[D_{\text{P}}\right]\left[E_{\text{P}}\right]{\left\{I_{\text{P}},R_{\text{P}},R_{\text{S}}\right\}}^{\text{T}},\end{array}

\begin{array}{}\text{(11b)}& {\left\{{\hat{\sigma }}_{xx},{\hat{\sigma }}_{yy},{\hat{\tau }}_{xy}\right\}}_{\text{P}}^{\text{T}}=\left[S_{\text{P}}\right]\left[E_{\text{P}}\right]{\left\{I_{\text{P}},R_{\text{P}},R_{\text{S}}\right\}}^{\text{T}},\end{array}

where

\begin{array}{}\text{(11c)}& \left[D_{\text{P}}\right]=\left[\begin{array}{ccc}-\text{i}{k}_{x\text{P}}& -\text{i}{k}_{x\text{P}}& -r_{\text{S}}\\ r_{\text{P}}& -r_{\text{P}}& \text{i}{k}_{x\text{P}}\end{array}\right],\end{array}

\begin{array}{}\text{(11d)}& \left[E_{\text{P}}\right]=\left[\begin{array}{ccc}\exp \left(r_{\text{P}}y\right)& 0& 0\\ 0& \exp \left(-r_{\text{P}}y\right)& 0\\ 0& 0& \exp \left(-r_{\text{S}}y\right)\end{array}\right],\end{array}

\begin{array}{}\text{(11e)}& \left[S_{\text{P}}\right]=\left[\begin{array}{ccc}-2\mu k_{x\text{P}}^2-\lambda k_{\text{P}}^2& -2\mu k_{x\text{P}}^2-\lambda k_{\text{P}}^2& 2\mu \mathrm{i}{k}_{x\text{P}}{r}_{\text{S}}\\ 2\mu r_{\text{P}}^2-\lambda k_{\text{P}}^2& 2\mu r_{\text{P}}^2-\lambda k_{\text{P}}^2& -2\mu \mathrm{i}{k}_{x\text{P}}{r}_{\text{S}}\\ -2\mu \mathrm{i}{k}_{x\text{P}}{r}_{\text{P}}& 2\mu \mathrm{i}{k}_{x\text{P}}{r}_{\text{P}}& \mu \left(r_{\text{S}}^2+k_{x\text{P}}^2\right)\end{array}\right].\end{array}

The free-field stress boundary conditions on the ground are:

\begin{array}{}\text{(12a)}& {\hat{\sigma }}_{yy}(y=0)=0,\end{array}

\begin{array}{}\text{(12b)}& {\hat{\tau }}_{xy}(y=0)=0.\end{array}

By substituting Eq. (12) into Eq. (11b), along with the definition for the wave velocity ratio:

(13)V=c_{\text{P}}/\phantom{c_{\text{P}}{c}_{\text{S}}}{c}_{\text{S}}=\sqrt{\left(\lambda +2\mu \right)/\phantom{\left(\lambda +2\mu \right)\mu }\mu },

one can obtain the amplitude ratios for the case with oblique incident P waves as:

\begin{array}{}\text{(14a)}& f_{\text{PP}}=\frac{R_{\text{P}}}{I_{\text{P}}}=-1-\frac{8k_{x\text{P}}^2{r}_{\text{P}}{r}_{\text{S}}}{\left(2k_{x\text{P}}^2-k_{\text{S}}^2\right)\left(2k_{x\text{P}}^2-V^2{k}_{\text{P}}^2\right)-4k_{x\text{P}}^2{r}_{\text{P}}{r}_{\text{S}}},\end{array}

\begin{array}{}\text{(14b)}& f_{\text{PS}}=\frac{R_{\text{S}}}{I_{\text{P}}}=\frac{4\text{i}{k}_{x\text{P}}{r}_{\text{P}}\left(2k_{x\text{P}}^2-V^2{k}_{\text{P}}^2\right)}{\left(2k_{x\text{P}}^2-k_{\text{S}}^2\right)\left(2k_{x\text{P}}^2-V^2{k}_{\text{P}}^2\right)-4k_{x\text{P}}^2{r}_{\text{P}}{r}_{\text{S}}},\end{array}

in which f_{\mathrm{P}\mathrm{P}} denotes the amplitude ratio of the reflected P waves to incident P waves and f_{\mathrm{P}\mathrm{S}} that of the reflected SV waves to the incident P waves.

Finally, by substituting the amplitude ratios into Eq. (11), the free-field responses of displacements and stresses in frequency and space domain can be solved by applying the inverse Fourier transformation with respect to k_{x\mathrm{P}}. The results in single Fourier transform (denoted with “-”) are

\begin{array}{}\text{(15a)}& {\left\{\overline{u},\overline{v}\right\}}_{\text{P}}^{\text{T}}=\left[D_{\text{P}}\right]\left[E_{\text{P}}\right]{\left\{1,f_{\text{PP}},f_{\text{PS}}\right\}}^{\text{T}}{I}_{\text{P}}\exp (-\text{i}{k}_{x\text{P}}x),\end{array}

\begin{array}{}\text{(15b)}& {\left\{{\overline{\sigma }}_{xx},{\overline{\sigma }}_{yy},{\overline{\tau }}_{xy}\right\}}_{\text{P}}^{\text{T}}=\left[S_{\text{P}}\right]\left[E_{\text{P}}\right]{\left\{1,f_{\text{PP}},f_{\text{PS}}\right\}}^{\text{T}}{I}_{\text{P}}\exp (-\text{i}{k}_{x\text{P}}x).\end{array}

The expressions given in Eqs. (15a) and (15b) represent the exact or closed-form solutions for the displacements and stresses, respectively, in frequency-space domain of the entire free field subjected to oblique P waves. Noteworthy is that, for the P waves with a fixed incident angle in both time and space domains, the wave number k_{x\mathrm{P}} remains constant in the process of wave propagation. Hence, one can multiply the responses in frequency and wavenumber domain by \mathrm{e} \mathrm{x}\mathrm{p}(- \mathrm{i}{k}_{x\mathrm{P}}x) to obtain the single inverse Fourier transform with respect to k_{x\mathrm{P}}.

3.2 Exact solution for incident SV waves

As revealed in Fig.3, for oblique SV-waves, the incident angle of the SV waves is denoted by {\theta }_{\mathrm{S}\mathrm{i}} and the reflected angle by {\theta }_{\mathrm{S}\mathrm{S}\mathrm{r}}, and the reflected angle of the P waves by {\theta }_{\mathrm{S}\mathrm{P}\mathrm{r}}. Similar to the case with P-wave incidence, combining the Snell’s law and the definitions of wave numbers for incident and reflected waves along x-axis (k_{x\mathrm{S}\mathrm{i}},k_{x\mathrm{S}\mathrm{P}\mathrm{r}},k_{x\mathrm{S}\mathrm{S}\mathrm{r}}), one can obtain the same relationship between the wave numbers as the case for oblique incident P waves, that is, k_{x\mathrm{S}\mathrm{i}} =k_{x\mathrm{S}\mathrm{P}\mathrm{r}}=k_{x\mathrm{S}\mathrm{S}\mathrm{r}}.

Let us define k_{x\mathrm{S}} \left(=\omega \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{S}\mathrm{i}}/c_{\mathrm{S}}\right) as the wave numbers of all waves along the x-axis for oblique SV-waves, i.e., with k_{x\mathrm{S}}=\omega \mathrm{s}\mathrm{i}\mathrm{n} {\theta }_{\mathrm{S}\mathrm{i}} /c_{\mathrm{S}}. The solutions of {\hat{\Phi }}_{\mathrm{S}} and {\hat{\Psi }}_{xy,\mathrm{S}}^{\mathrm{S}\mathrm{V}} in Eq. (9) are

\begin{array}{}\text{(16a)}& {\hat{\Phi }}_{\text{S}}=R_{\text{P}}\exp \left(-r_{\text{P}}y\right),\end{array}

\begin{array}{}\text{(16b)}& {\hat{\Psi }}_{xy\text{,S}}^{\text{SV}}=I_{\text{S}}\exp \left(r_{\text{S}}y\right)+R_{\text{S}}\exp \left(-r_{\text{S}}y\right),\end{array}

where r_{\mathrm{P}}^2=k_{x\mathrm{S}}^2-k_{\mathrm{P}}^2, r_{\mathrm{S}}^2= k_{x\mathrm{S}}^2- k_{\mathrm{S}}^2, I_{\mathrm{S}} represents the amplitude of the incident SV waves and R_{\mathrm{S}} that of the reflected ones, and R_{\mathrm{P}} is the amplitude of the reflected P waves. The direction of propagation of each wave component is consistent with that of the incident P waves.

Substituting the potential function in Eq. (16) into Eqs. (6) and (7) yields the displacements and stresses in terms of R_{\mathrm{P}}, I_{\mathrm{S}}, and R_{\mathrm{S}} in transformed form as:

\begin{array}{}\text{(17a)}& {\left\{\hat{u},\hat{v}\right\}}_{\text{S}}^{\text{T}}=\left[D_{\text{S}}\right]\left[E_{\text{S}}\right]{\left\{R_{\text{P}},I_{\text{S}},R_{\text{S}}\right\}}^{\text{T}},\end{array}

\begin{array}{}\text{(17b)}& {\left\{{\hat{\sigma }}_{xx},{\hat{\sigma }}_{yy},{\hat{\tau }}_{xy}\right\}}_{\text{S}}^{\text{T}}=\left[S_{\text{S}}\right]\left[E_{\text{S}}\right]{\left\{R_{\text{P}},I_{\text{S}},R_{\text{S}}\right\}}^{\text{T}},\end{array}

where

\begin{array}{}\text{(17c)}& \left[D_{\text{S}}\right]=\left[\begin{array}{ccc}-\text{i}{k}_{x\text{S}}& r_{\text{S}}& -r_{\text{S}}\\ -r_{\text{P}}& \text{i}{k}_{x\text{S}}& \text{i}{k}_{x\text{S}}\end{array}\right],\end{array}

\begin{array}{}\text{(17d)}& \left[E_{\text{S}}\right]=\left[\begin{array}{ccc}\exp \left(-r_{\text{P}}y\right)& 0& 0\\ 0& \exp \left(r_{\text{S}}y\right)& 0\\ 0& 0& \exp \left(-r_{\text{S}}y\right)\end{array}\right],\end{array}

\begin{array}{}\text{(17e)}& \left[S_{\text{S}}\right]=\left[\begin{array}{ccc}-2\mu k_{x\text{S}}^2-\lambda k_{\text{P}}^2& -2\mu \text{i}{k}_{x\text{S}}{r}_{\text{S}}& 2\mu \text{i}{k}_{x\text{S}}{r}_{\text{S}}\\ 2\mu r_{\text{P}}^2-\lambda k_{\text{P}}^2& 2\mu \text{i}{k}_{x\text{S}}{r}_{\text{S}}& -2\mu \text{i}{k}_{x\text{S}}{r}_{\text{S}}\\ 2\mu \text{i}{k}_{x\text{S}}{r}_{\text{P}}& \mu \left(r_{\text{S}}^2+k_{x\text{S}}^2\right)& \mu \left(r_{\text{S}}^2+k_{x\text{S}}^2\right)\end{array}\right].\end{array}

With Eq. (12) substituted into Eq. (17b), the amplitude ratios caused by oblique incident SV waves can be obtained as:

\begin{array}{}\text{(18a)}& f_{\text{SP}}=\frac{R_{\text{P}}}{I_{\text{S}}}=\frac{-4\text{i}{k}_{x\text{S}}{r}_{\text{S}}\left(2k_{x\text{S}}^2-k_{\text{S}}^2\right)}{\left(2k_{x\text{S}}^2-k_{\text{S}}^2\right)\left(2k_{x\text{S}}^2-V^2{k}_{\text{P}}^2\right)-4k_{x\text{S}}^2{r}_{\text{P}}{r}_{\text{S}}},\end{array}

\begin{array}{}\text{(18b)}& f_{\text{SS}}=\frac{R_{\text{S}}}{I_{\text{S}}}=-1-\frac{8k_{x\text{S}}^2{r}_{\text{P}}{r}_{\text{S}}}{\left(2k_{x\text{S}}^2-k_{\text{S}}^2\right)\left(2k_{x\text{S}}^2-V^2{k}_{\text{P}}^2\right)-4k_{x\text{S}}^2{r}_{\text{P}}{r}_{\text{S}}},\end{array}

in which f_{\mathrm{S}\mathrm{P}} denotes the amplitude ratio of the reflected P waves to the incident SV waves and f_{\mathrm{S}\mathrm{S}} the amplitude ratio of the reflected SV waves to the incident SV waves.

Finally, substituting the amplitudes into Eq. (17), together with the inverse Fourier transformation with respect to k_{x\mathrm{S}}, yields the expression of the displacements and stresses (in single Fourier transform) of the free field in frequency and space domain:

\begin{array}{}\text{(19b)}& {\left\{\overline{u},\overline{v}\right\}}_{\text{S}}^{\text{T}}=\left[D_{\text{S}}\right]\left[E_{\text{S}}\right]{\left\{f_{\text{SP}},1,f_{\text{SS}}\right\}}^{\text{T}}{I}_{\text{S}}\exp (-\text{i}{k}_{x\text{S}}x),\end{array}

\begin{array}{}\text{(19b)}& {\left\{{\overline{\sigma }}_{xx},{\overline{\sigma }}_{yy},{\overline{\tau }}_{xy}\right\}}_{\text{S}}^{\text{T}}=\left[S_{\text{S}}\right]\left[E_{\text{S}}\right]{\left\{f_{\text{SP}},1,f_{\text{SS}}\right\}}^{\text{T}}{I}_{\text{S}}\exp (-\text{i}{k}_{x\text{S}}x).\end{array}

Similar to the case with P-wave incidence, the expressions given in Eqs. (19a) and (19b) represent the exact or closed-form solutions for the displacements and stresses, respectively, in frequency-space domain of the entire free field subjected to oblique incident SV waves.

From Fig.3, it was observed that the reflected angle of P waves under any condition is greater than that of SV waves for oblique SV-wave incidence. What is more, a critical angle {\theta }_{\mathrm{s}\mathrm{i}}^{\mathrm{c}\mathrm{r}} exists for the incident SV waves, above which the reflected P waves will disappear. Under such a condition, the reflected P waves caused by incident SV waves begin to transform from body waves to surface waves. For instance, by letting \mathrm{s}\mathrm{i}\mathrm{n} {\theta }_{\mathrm{S}\mathrm{P}\mathrm{r}}=1, along with the wave velocity ratio defined as V=c_{\mathrm{P}}/c_{\mathrm{S}}, the critical angle of the incident SV waves can be obtained from Fig.3 as:

(20){\theta }_{\text{Si}}^{\text{cr}}=\arcsin \frac{1}{V}.

In the previous works [2,3], the potential functions and dynamic responses were both given in time domain, which requires the relevant equations for SV waves to be divided into two parts for incident angles below or above the critical one, but not in a unified format. This makes the implementation of the ground motion input in numerical modeling more difficult. Such a division is unnecessary when using the frequency domain solution method, as will be discussed below concerning the applicability of the response of free field for the SV waves in Eq. (19).

Firstly, according to the relationship between the incident and reflected angles in Fig.3, for {\theta }_{\mathrm{S}\mathrm{i}} >{\theta }_{\mathrm{S}\mathrm{i}}^{\mathrm{c}\mathrm{r}}, it is required that \mathrm{s}\mathrm{i}\mathrm{n} {\theta }_{\mathrm{S}\mathrm{P}\mathrm{r}}=1, where {\theta }_{\mathrm{S}\mathrm{P}\mathrm{r}} does not really exist in the geometric space, but should be interpreted as a complex angle following Lee and Karl [5]. Then, based on the wave numbers definition k_{x\mathrm{S}}=k_{x\mathrm{S}\mathrm{P}\mathrm{r}}=\omega \mathrm{s}\mathrm{i}\mathrm{n} {\theta }_{\mathrm{S}\mathrm{P}\mathrm{r}} /c_{\mathrm{P}}, one also finds that k_{x\mathrm{S}}>k_{\mathrm{P}}, since \mathrm{s}\mathrm{i}\mathrm{n} {\theta }_{\mathrm{S}\mathrm{P}r}=1, and accordingly r_{\mathrm{P}}>0 is valid, as indicated below:

(21)r_{\text{P}} =\sqrt{k_{x\text{S}}^2-k_{\text{P}}^2}=\sqrt{{\left(\frac{\omega }{c_{\mathrm{P}}}\right)}^2\left({\sin }^2{\theta }_{\text{SPr}}-1\right)}>0.

Subsequently, for r_{\mathrm{P}}>0, the power \left({-r}_{\mathrm{P}}y\right) in the exponent \mathrm{e}\mathrm{x}\mathrm{p} \left({-r}_{\mathrm{P}}y\right) of the potential function of the reflected P waves in Eq. (16a) becomes positive real for \left(-y\right) (where y is positive downward), which satisfies the condition of surface waves, in that the wave decays rapidly with increasing depth as noted by Rayleigh [58]. Thus, the reflected P waves appears as surface waves, rather than body waves, traveling parallel to the surface with the amplitude decaying along the depth direction (y-axis), as shown in Fig.4. Such a result is consistent with Lee and Karl [5] using the Fourier-Bessel series in the polar coordinates. Hence, the potential function in Eq. (16a) of the reflected P waves R_{\mathrm{P}}\mathrm{e}\mathrm{x}\mathrm{p} \left({-r}_{\mathrm{P}}y\right) in frequency and wavenumber domain can be used for oblique incident SV waves. Unlike Lee and Karl’s [5] formulation in the polar coordinates, the dynamic responses obtained from the potential (see Eq. (19)) in the Cartesian coordinates in frequency and space domain are valid for all SV waves with arbitrary incident angles, including the critical angle. This has made the ground motion input rather easy and straightforward, as will be described in Sections 4 and 5 to follow.

The amplitudes \left(I_{\mathrm{P}},I_{\mathrm{S}}\right) of incident waves involved in Eqs. (15) and (19) will be determined by additional boundary conditions in Subsection 4.1 for the acceleration responses. Moreover, the amplitude ratios in Eqs. (14) and (18) under oblique P- and SV-waves deduced from the general solutions of the potential functions in frequency domain in this paper are expressed only in terms of the wave velocity ratio and wave numbers, irrelevant of the reflected angles as was the case with Ewing and Graff’s work [3]. Meanwhile, the solutions for the reflected angles do not directly apply to the incidence SV waves beyond the critical angle.

In short, the present formulation for the dynamic responses of the free field in closed form not only simplifies the procedure of derivation, but also provides a rational basis for inclusion of oblique incident P and SV waves in the soil dynamic analysis. Compared with previous works, the closed-form solutions derived for 2D oblique P- and SV-waves are neater in frequency domain, especially applicable for SV waves with over-critical incident angles. With the solution presented above, the frequency-based equivalent seismic forces for use in the finite/infinite element analysis can be computed in a straightforward way, as will be shown in the following section. It should be added that the above procedure for dealing with oblique incident SV-waves can be applied to layered soils as well, which is reserved for a future work.

4 Equivalent seismic forces under oblique incident seismic waves

As schematically illustrated in Fig.5, two regions are considered for the half-space, i.e., near field (containing the part of soil of interest) and the far field (extending to infinity). For convenience, the near field is taken as a rectangular region, in the numerical analysis to follow. The symbols b_L{,b}_R{,b}_B are used to denote the interfacing boundaries on the left, right and bottom sides of the near field, respectively. In this paper, the near field will be discretized and modeled by finite elements and the far field by infinite elements, which will be briefed later on. For numerical analysis, one needs to determine the nodal load vector or equivalent seismic forces \mathit{F} exerted by oblique incident seismic body waves (from the far field) onto the interfacing boundary \mathit{b}, based on the available recorded data of the earthquake on the free surface, as will be presented below.

4.1 Displacements and stresses on near-field boundary

Before the equivalent seismic forces can be calculated, the displacements and stresses of the near-field boundary caused by unit ground acceleration should first be determined. For oblique incident seismic waves, a certain angle exists between the wave front and the ground. As shown in Fig.6, it is assumed that each point on the near-field boundary is subjected to both incident and reflected waves. Accordingly, both the incident and reflected waves should be included in calculation of the responses for the near-field boundary, in consistence with the theoretical formula for the free-field response presented in Subsections 3.1 and 3.2.

The earthquake to hit from the far field is assumed to have been recorded (and made available) at the origin on the free-field surface. In order to transform the recorded seismic data into the response of the near-field boundary, one needs to determine for each frequency \omega how a unit horizontal acceleration of incident SV waves or a unit vertical acceleration of incident P waves on the ground is transformed to the near-field boundary. Namely, consider the unit accelerations acting at the origin of the free surface:

\begin{array}{}\text{(22a)}& -{\omega }^2\hat{u}(y=0)=1,\end{array}

\begin{array}{}\text{(22b)}& -{\omega }^2\hat{v}(y=0)=1.\end{array}

Accordingly, the displacements (horizontal or vertical) on the ground in frequency and wavenumber domain for incident SV or P waves, respectively, are

\begin{array}{}\text{(23a)}& \hat{u}(y=0)=-1/\phantom{-1{\omega }^2}{\omega }^2,\end{array}

\begin{array}{}\text{(22b)}& \hat{v}(y=0)=-1/\phantom{-1{\omega }^2}{\omega }^2.\end{array}

For incident P waves, one can substitute the displacement \hat{v} into Eq. (11a), and for incident SV waves, the displacement \hat{u} into Eq. (17a), along with the amplitude ratios in Eqs. (14) and (18), respectively, to obtain the amplitudes of the two incident waves in frequency-wavenumber domain as:

\begin{array}{}\text{(24a)}& I_{\text{P}}=-{\left[{\omega }^2\left(r_{\text{P}}-r_{\text{P}}{f}_{\text{PP}}+\text{i}{k}_{x\text{P}}{f}_{\text{PS}}\right)\right]}^{-1},\end{array}

\begin{array}{}\text{(24b)}& I_{\text{S}}=-{\left[{\omega }^2\left(-\text{i}{k}_{x\text{S}}{f}_{\text{SP}}+r_{\text{S}}-r_{\text{S}}{f}_{\text{SS}}\right)\right]}^{-1}.\end{array}

Subsequently, for each point \left(x,y\right) on the near-field boundary in Fig.5, one can substitute the amplitudes of the incident waves in Eq. (24) into the general solutions in Eqs. (15) and (19) for P- and SV-waves, respectively, to obtain the free-field displacements {\mathit{u}}_{\mathrm{b}}^{\mathrm{f}} and stresses {\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} in frequency and space domain for each unit acceleration of frequency \omega. Such a procedure can be looped over for all the spectral frequencies covered by each seismic record, to generate the corresponding displacements {\mathit{u}}_{\mathrm{b}}^{\mathrm{f}} and stresses {\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} on the near-field boundary. Based on this, the equivalent seismic forces can be calculated, as will be presented in the section to follow.

4.2 Equivalent seismic forces on near-field boundary

The equivalent seismic forces imposed on the near-field boundary nodes will be calculated from the displacements {\mathit{u}}_{\mathrm{b}}^{\mathrm{f}} and stresses {\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} (mentioned in Subsection 4.1) for a given free field seismic record. For illustration, the acceleration history of the Taft Earthquake recorded at No.1095 station in Lincoln University in 1952 with the sampling interval \mathrm{\Delta }t= 0.02 s will be adopted. The horizontal acceleration with peak ground acceleration (PGA) of 155.7 gal in Fig.7(a) is chosen as the one for incident SV waves, and the vertical acceleration with peak ground acceleration (PGA) of 104.9 gal in Fig.7(b) for incident P waves. Subsequently, by performing the fast Fourier transform (FFT), the acceleration frequency spectra for both directions were plotted in Fig.8, where the horizontal axis denotes the circular frequency in Hz (= 2π rad/s). In general, to ensure proper resolution in frequency domain, the total number N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} of frequency increments of the FFT should be equal to the number of sampling points for the acceleration time histories. According to Nyquist’s theorem, since the frequency range from 0 to 1/\left(2\mathrm{\Delta }t\right) is the one valid for sampling, the total number of frequency increments selected for generating the spectrum is half of the total number N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} mentioned above, denoted by N_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}}, namely, the frequency increment used for generating the spectrum in the FFT is: \mathrm{\Delta }\omega =\text{π}/\left[\mathrm{\Delta }t\left(N_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}}-1\right)\right] rad/s.

As can be seen from Fig.8, the earthquake excitation is concentrated mainly in the low frequency part, namely, it become negligible for frequencies over 25 Hz. Thus, only frequencies up to 25 Hz will be included for the ground motion input.

After obtaining the acceleration spectra in Fig.8, the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} for each frequency \omega can be related to the boundary displacements {\mathit{u}}_{\mathrm{b}}^{\mathrm{f}} and stresses {\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} following Zhao and Valliappan [47]. Accordingly, the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} on the left boundary b_L are

\begin{array}{}\text{(25a)}& {\mathit{F}}_{\text{eq}}={\mathit{S}}_{\text{rb}}{\mathit{u}}_{\text{b}}^{\text{f}}-{\mathit{A}\mathit{\sigma }}_{\text{b}}^{\text{f}},\end{array}

and those on the right and bottom boundaries \left(b_{\mathrm{R}},b_{\mathrm{B}}\right) are

\begin{array}{}\text{(25b)}& {\mathit{F}}_{\text{eq}}={\mathit{S}}_{\text{rb}}{\mathit{u}}_{\text{b}}^{\text{f}}+{\mathit{A}\mathit{\sigma }}_{\text{b}}^{\text{f}},\end{array}

where {\mathit{S}}_{\mathrm{r}\mathrm{b}} is the dynamic impedance matrix of the infinite elements used to represent the actual effect of the far field on the boundary (to be elaborated in Section 5), and A the area matrix for transforming the surface tractions into nodal forces. Since the nodal stresses generated by seismic waves on the left boundary are in the opposite direction to those on the right and bottom with respect to the x-axis, the surface tractions \mathit{A}{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{f}} in Eqs. (25a) and (25b) differ by a negative sign. In the numerical analysis to follow, the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} given in Eqs. (25a) and (25b) will be imposed as the external loads on the near-field boundary.

5 Seismic analysis by finite/infinite element method

In this section, the procedure for simulating the half-space under oblique seismic waves by the FIEM will be given. In this case, the oblique incident seismic waves hitting the near field are asymmetric with respect to the vertical central line. For convenience, the origin O is set on the left-top corner of the finite/infinite element mesh, see Fig.9. The finite elements will be adopted to simulate the near field and infinite elements the far field. In this study, the infinite elements have been arranged to have oblique infinite boundaries consistent with the wave propagation directions for treating the oblique incident seismic waves, as shown in Fig.9. The width and depth of the near field (i.e., finite element mesh) are taken as 60 m and 30 m along the x- and y -axes, respectively. Here, the observation points (OP-1−OP-6) in Fig.9 are the ones to be referred to in the numerical analyses.

5.1 Equation of motion by finite/infinite element method

The elements adopted are the quadratic 8-node (Q8) finite element and the degenerated 5-node (Q5) infinite element [48]. For brevity, only the shape functions for the displacements of the Q5 infinite elements (see Fig.10) will be briefed as:

\begin{array}{}\text{(26a)}& N_1=\frac{\eta \left(\eta -1\right)}{2}P\left(\xi \right),\end{array}

\begin{array}{}\text{(26b)}& N_2=-\left(\eta -1\right)\left(\eta +1\right)P\left(\xi \right),\end{array}

\begin{array}{}\text{(26c)}& N_3=\frac{\eta \left(\eta +1\right)}{2}P\left(\xi \right),\end{array}

(27)P\left(\xi \right)={\text{e}}^{-{\alpha }_{\text{l}}\xi }{\text{e}}^{-\text{i}{k}_{\text{l}}\xi },

where P\left(\xi \right) represents the wave propagation function in the local coordinates, {\alpha }_{\mathrm{l}} denotes the amplitude decay factor with the decaying function {\mathrm{e}}^{-{\alpha }_{\mathrm{l}}\xi } representing the radiation damping induced by waves propagating to infinity, and k_{\mathrm{l}} the wave number with {\mathrm{e}}^{-{\mathrm{i}k}_{\mathrm{l}}\xi } representing the effect of wave oscillation. The function P\left(\xi \right) in Eq. (27) satisfies the Sommerfeld radiation condition, in that it can represent the decay of outward traveling waves to zero at infinity [59]. Selection of the parameters {\alpha }_{\mathrm{l}} and k_{\mathrm{l}} are crucial to the employment of infinite elements, which will be briefed in Subsection 5.2. The transformation of the wave propagation function from the local coordinates to the global ones is available in Ref. [48], which will not be recapitulated herein.

Consider that the external forces and displacements of the nodes are of the harmonic form. According to the standard procedure of the FEM, the motion equations for the half-space can be established for each frequency \omega as:

(28)\left(\mathit{K}-{\omega }^2\mathit{M}\right)\mathit{U}= \mathit{F},

where K and M denote the stiffness and mass matrices, respectively [48], for the discretized half-space shown in Fig.5, \mathit{U} is the nodal displacements vector, and F the nodal load vector in the global coordinates. The matrix \left(\mathit{K}-{\omega }^2\mathit{M}\right) is also known as the impedance matrix \mathit{S}\left(\omega \right).

In soil dynamic analysis, one may divide the impedance matrix \mathit{S}\left(\omega \right) into two parts related to the near field (n) and boundary (b) and recast Eq. (28) as follows:

(29)\left[\begin{array}{cc}{\mathit{S}}_{\text{nn}}\left(\omega \right)& {\mathit{S}}_{\text{nb}}\left(\omega \right)\\ {\mathit{S}}_{\text{bn}}\left(\omega \right)& {\mathit{S}}_{\text{bb}}\left(\omega \right)+{\mathit{S}}_{\text{rb}}\left(\omega \right)\end{array}\right]\left[\begin{array}{c}{\mathit{U}}_{\text{n}}\left(\omega \right)\\ {\mathit{U}}_{\text{b}}\left(\omega \right)\end{array}\right] =\left[\begin{array}{c}0\\ {\mathit{F}}_{\text{eq}}(\omega )\end{array}\right],

where {\mathit{F}}_{\mathrm{e}\mathrm{q}} \left(\omega \right) is exactly the equivalent seismic forces obtained in Section 4 for a unit ground acceleration, {\mathit{S}}_{\mathrm{n}\mathrm{n}} \left(\omega \right),{\mathit{S}}_{\mathrm{n}\mathrm{b}} \left(\omega \right),{\mathit{S}}_{\mathrm{b}\mathrm{n}} \left(\omega \right),{\mathit{S}}_{\mathrm{b}\mathrm{b}} \left(\omega \right) represent the near-field impedance matrices and interfacing boundary, and {\mathit{U}}_{\mathrm{n}}\left(\omega \right),{\mathit{U}}_{\mathrm{b}}\left(\omega \right) are the near-field displacements and boundary, respectively, in frequency and space domain.

It is worth noting that infinite elements serve to absorb the scattered waves generated by any excitation imposed on the near field, including the equivalent seismic forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} imposed on the near-field boundary. The capability of the infinite elements in dealing with the geometry radiation of the far field for waves propagating toward infinity has been fully verified with the key parameters specified in Ref. [48], which will be modified for the problem below.

5.2 Amplitude decay factors {\alpha }_{\mathrm{l}} and wave numbers k_{\mathrm{l}}

The key parameters of the numerical elements need to be discussed to ensure the accuracy of calculation. First, let us talk about the amplitude decay factor {\alpha }_{\mathrm{l}} used for the infinite elements in Eq. (27). For seismic waves with incident angles in the range , both P and SV waves will propagate in near field, and Rayleigh (R) waves will also occur on the ground due to the interference of body waves. The present problem for seismic waves with incident angles oblique to the half-space with symmetrical near-field region is similar to the wave propagation problem led by a line load acting at the center point C of the near-field ground in Fig.11, in which Q denotes the distance between the center point C of the ground and the near-field boundary. According to Ref [48], for the regions where the body waves (i.e., P and SV waves) are dominant, the amplitude decay factor {\alpha }_{\mathrm{l}} can be taken as 1 /2 Q, and for the regions where R waves are dominant, {\alpha }_{\mathrm{l}} is taken as zero for the 2D problem considered. The above treatment remains valid for the case for SV waves with incident angles over the critical one.

Next, let us discuss the selection of wave numbers used for the infinite elements in Eq. (27). Two different partitioning schemes need to be considered. For the case with P waves and SV waves with incident angles below the critical angle, the R- and SV-wave numbers, k_{\mathrm{R}} and k_{\mathrm{S}}, respectively, are used for the infinite elements in the region near and below the ground on two sides, and P-wave number k_{\mathrm{P}} for the region at deep bottom, as shown in Fig.11(a) [48]. But, for the case with SV waves with incident angles over the critical angle, the SV-wave number k_{\mathrm{S}} is used for both the region below the ground on two sides and at deep bottom, due to absence of P waves. The accuracy of the above selection of wave numbers will be confirmed in the numerical verification.

5.3 Procedure of analysis for oblique incident seismic waves

The seismic analysis procedure for oblique incident body waves, especially the transformation of the solutions from frequency to time domain by the inverse fast Fourier transform (IFFT), will be outlined in this section.

Firstly, the range of frequencies considered, i.e., 25 Hz (see Fig.8), for the earthquake is divided into a number of increments for frequency domain calculation. Here, an increment of \mathrm{\Delta }\omega = 0.1 Hz (= 0.6283 rad/s) has been proved to be sufficient for the accuracy of calculation [54], and will be adopted throughout. It should be added that the frequency increment used for generating the earthquake spectra in Fig.8 using the FFT is approximately 0.015 Hz.

Then, for each frequency \omega within the range, the equivalent boundary forces {\mathit{F}}_{\mathrm{e}\mathrm{q}} \left(\omega \right) can be computed by Eq. (25). The displacements \mathit{U}\left(\omega \right) of the near field and boundary can be solved from Eq. (29), and multiplied by -{\omega }^2 to yield the acceleration responses {\mathit{U}}^{″}\left(\omega \right). By Hooke’s law, one can also calculate the stresses \mathit{\sigma }\left(\omega \right) for all nodes of the near field in frequency and space domain. The set of component responses for the displacements, accelerations and stresses are those obtained for the frequency \omega. The procedure can be repeated for the next frequency \omega +\mathrm{\Delta }\omega. It should be noted that each set of component responses has a corresponding spectral amplitude in the earthquake spectra of Fig.8.

Next, by multiplying the above set of component responses for each frequency \omega by the corresponding spectral amplitude of the earthquake spectra in Fig.8, one can obtain the set of spectral responses contributed by each frequency component \omega of the earthquake for the horizontal (SV) or vertical (P) motion in Fig.7.

There are two integrals for the IFFT in Eq. (8b). For the integral with respect to the wavenumber (k_{x\mathrm{P}} or k_{x\mathrm{S}}), it was done in calculating the theoretical solutions of the displacement and stress responses of the near-field boundary (see Eqs. (15) and (19)) used in the equivalent seismic forces computation in Eq. (25). Hence, after obtaining the near-field responses in frequency and space domain (see Eq. (29)), only the integral with respect to frequency \omega in Eq. (8b) needs to be implemented. In practice, this is actually performed by replacing the integral by the summation of all the spectral responses for all frequencies falling within the range of consideration for each earthquake [54]. In other words, by looping over the frequency \omega from 0 to 25 Hz with the increment \mathrm{\Delta }\omega (rad/s) for the IFFT, one can sum up all the spectral responses to obtain the total time-history responses of the whole region of interest, including near field and boundary, for which the time increment \mathrm{\Delta }t is identical to the one adopted in creating the time-hisotry records for the earthquake. The numerical procedure for seismic analysis is implemented by the MATLAB software herein.

In this paper, the analysis of seismic response based on the frequency-domain superposition method saves the tedious computation for the step-by-step analysis in time domain. Using the present method, the calculation efficiency is also improved for analysis of the dynamic response of the site and underground structures considered herein under the action of oblique incident earthquakes.

6 Verification of seismic analysis procedure

6.1 Verification of frequency domain solution with Wolf’s

In the numerical verification, the finite element size L is set to be 1 m × 1 m, and the mesh range of Q = 30 m is selected for the model in Fig.9. To compare the present solutions with Wolf’s [4], the physical properties of the soil adopted are: elastic modulus E_{\mathrm{s}}= 100 MPa, Poisson’s ratio \nu = 0.33, the density {\rho }_{\mathrm{s}}= 2000 kg/m³, for both P- and SV-wave incidence. Hysteresis damping was not considered in Wolf’s work [4]. For comparison, it will be neglected here too, but will be included in subsequent numerical analyses. According to Eq. (20), the critical angle of the SV waves is 30°.

For the earthquake spectra shown in Fig.8, the range of frequency considered is up to 25 Hz. By applying a unit vertical or horizontal acceleration of incident P or SV waves at the origin (see Fig.9) for each frequency \omega (see Subsection 4.1), the normalized horizontal and vertical ground displacement amplitudes \left|\overline{u}/U_{\mathrm{P}}\right|, \left|\overline{v}/U_{\mathrm{P}}\right|, \left|\overline{u}/U_{\mathrm{S}}\right|, and \left|\overline{v}/U_{\mathrm{S}}\right|, computed by present method, along with those in Wolf’s were plotted in Fig.12 and Fig.13, for both P- and SV-waves, respectively. Here, U_{\mathrm{P}} and U_{\mathrm{S}} are the total ground displacement amplitudes caused by upward P and SV waves, respectively.

As can be observed from Fig.12 and Fig.13, all the numerical solutions agree excellently with the analytical solutions of Wolf for incident angles from 0° to 90°. Of interest is the drastically changing behavior of the soil around the critical angle of 30° in Fig.13(a) and 13(b). It is confirmed that the infinite elements adopted for the near-field boundaries \left(b_{\mathrm{L}}{,b}_{\mathrm{R}}{,b}_{\mathrm{B}}\right) can absorb adequately the inhomogeneous surface waves generated by the SV waves incident waves around the critical angle of 30°. Besides, it is inferred that the element size and mesh range selected herein (L = 1 m, Q = 30 m) are adequate, as the requirement of L\le {\lambda }_{\mathrm{S}}/5 has been met [55]. Consequently, the same element size and mesh range will be used in the subsequent simulations for the sake of accuracy.

6.2 Verification of frequency domain solution with Luco and de Barros’

For the half-space containing a cylindrical cavity (see Fig.14), the present solutions will be compared with de Barros and Luco’s [11] to verify the capability of infinite elements in absorbing the scattering waves generated by both P- and SV-waves. Following Ref. [11], the embedment depth is H = 10 m, and radius of cylindrical cavity is a = 5 m. The other properties of the soil are: Poisson’s ratio \nu = 0.33, elastic modulus E_{\mathrm{s}}= 100 MPa, density {\rho }_s= 2000 kg/m³, hysteresis damping \beta = 0.01.

For comparison with de Barros and Luco’s solutions [11], the dimensionless frequency \eta =2a/{\lambda }_{\mathrm{S}} is taken to be 0.5. The unit vertical or horizontal acceleration of incident P or SV waves, both with incident angle 45°, is applied at the point (30,0) in Fig.14. Using the same mesh range Q and finite element size L in Subsection 6.1, the normalized horizontal and vertical ground displacement amplitudes \left|\overline{u}/U_{\mathrm{P}}\right|, \left|\overline{v}/U_{\mathrm{P}}\right|, \left|\overline{u}/U_{\mathrm{S}}\right|, and \left|\overline{v}/U_{\mathrm{S}}\right| computed by present method, along with those of de Barros and Luco’s [11], were plotted in Fig.15 and Fig.16, for P- and SV-waves, respectively. As can be seen, the present solutions agree excellently with those of de Barros and Luco’ for P- and SV-waves. This has verified the capability of the infinite elements in absorbing the scattered waves generated by oblique P- and SV-waves. On the whole, the use of finite/infinite elements makes it easy not only to simulate the far field with high accuracy, but also the half-space containing underground structures in a simple way [53–55], which are the advantages over the other methods.

6.3 Verification of time domain solution

Unless otherwise specified, the hysteretic damping ratio is taken as \beta = 0.02, Poisson’s ratio as \nu = 0.25, the density as {\rho }_{\mathrm{s}}= 2000 kg/m³, and the elastic modulus as E_{\mathrm{s}}= 100 MPa for the following numerical analysis. Aimed at verifying the accuracy of the treatment for transformation of the solutions from the frequency to time domain for oblique incident seismic waves, as presented in the preceding section, the incident angles {\theta }_{\mathrm{P}\mathrm{i}} and {\theta }_{\mathrm{S}\mathrm{i}} (see Fig.2 and Fig.3) are set to 50°. Point OP-1 (on the surface) in Fig.9 was taken as the observation station, which was also the point where the Taft Earthquake was originally applied for analysis.

Both the inversely computed and original time-histories of the Taft Earthquake (see Fig.6 for original) at point OP-1 for the SV- and P-waves were, respectively, plotted in Fig.17(a) and Fig.17(b). The excellent agreement between the two sets of earthquake data is an indication that the procedure adopted for the IFFT for computing the time-history soil responses is accurate.

7 Parametric study of half-space

7.1 Frequency domain responses of half-space for P-wave incidence

Using the same numerical model in Fig.9 and the data previously given, the normalized horizontal and vertical displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) obtained for the free surface were plotted against frequency in Fig.18, and those at specified depths along the vertical central line of the soil model in Fig.19. For oblique incident P waves, one observes from Fig.18 that the ground attenuations \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) are more obvious along the x-axis for increasing frequency, implying that the ground responses far from the epicenter are mainly controlled by lower frequency seismic waves. However, from Fig.19, it is observed that the underground displacements \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) fluctuate with increasing frequency, revealing some local maxima and minima due to the soil resonance effect. Moreover, the number of resonant frequencies increases with respect to increasing depth, which agrees with Wolf’s conclusions [4].

With regard to the effect of incident angle at different depths, the analysis was focused on the displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) of the observation points along the vertical central line for incident angles {\theta }_{\mathrm{P}\mathrm{i}} varying from 0° to 90° under f = 25 Hz. As can be observed from Fig.20, there exist several local maxima and minima with the underground displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) for increasing incident angle at different depths, and the fluctuation appears to be more obvious at larger depth, different from the ground responses in Fig.12. It is concluded that the displacement responses at larger depth are more sensitive to the change in incident angle.

To further investigate the spatially irregular responses of the half-space, the horizontal and vertical displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) along x- and y-axes fixed at certain frequency and incident angle, i.e., f = 25 Hz, {\theta }_{\mathrm{P}\mathrm{i}} = 50°, were plotted in Fig.21. Here, to show the responses along the depth direction more clearly, the horizontal axis in Fig.21(b) represents the normalized horizontal or vertical displacement amplitudes, which will also be adopted in the next sections. From Fig.21(a), it is observed that both the horizontal and vertical displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) monotonously decrease along the horizontal direction of propagation for incident waves. As for the horizontal and vertical displacement amplitude \left(\left|\overline{u}/U_{\mathrm{P}}\right|,\left|\overline{v}/U_{\mathrm{P}}\right|\right) shown in Fig.21(b), they appear to be fluctuating synchronously along the depth direction with little attenuation and the displacement amplitudes at some certain depths are greater than those on free surface.

7.2 Frequency domain responses of half-space for SV-wave incidence

For oblique SV-wave incidence, the critical angle calculated of homogeneous soil is 35.26° by Eq. (20). With the incident angle {\theta }_{Si} = 50°, the normalized horizontal and vertical displacements \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) on free surface and at specified depths on vertical central line in frequency domain were, respectively, plotted in Fig.22 and Fig.23. As can be seen, the trends of displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) on free surface and at specified depths for increasing frequency are similar for incident P waves. One distinct difference is that the vertical displacement amplitudes \left(\left|\overline{v}/U_{\mathrm{S}}\right|\right) are significantly higher than the horizontal ones \left(\left|\overline{u}/U_{\mathrm{S}}\right|\right) especially on the ground for SV waves with over-critical incident angles, which is consistent with the phenomenon revealed in Fig.13.

Subsequently, the underground displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) on the vertical central line of the soil model for incident SV waves with incident angle {\theta }_{Si} varying from 0° to 90° for f = 25 Hz plotted in Fig.24. Clearly, the same trend in responses as the one for incident P waves can be observed. But the fluctuation of displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) with angle variation are significantly greater than that for incident P waves, indicating that the responses of displacement amplitudes have higher sensitivity to the change of incident angle under incident SV waves. Moreover, abrupt changes in both horizontal and vertical displacement amplitudes occur near the critical angle (see red lines), but they are not as obvious as the ground response in Fig.13.

In addition, for {\theta }_{Si}= 50° and f = 25 Hz, the horizontal and vertical displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) along the free surface and the vertical central line of the soil model were, respectively, plotted in Fig.25(a) and 25(b). As can be seen, except for the same trend as for P-wave incidence, the deviation between the vertical and horizontal displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) under SV waves with over-critical incident angles is significantly larger than that for incident P waves in Fig.21(a). Moreover, with increasing depth, the overall trend of horizontal and vertical displacement amplitudes \left(\left|\overline{u}/U_{\mathrm{S}}\right|,\left|\overline{v}/U_{\mathrm{S}}\right|\right) in Fig.24(b) differs by half a period of fluctuation, also different from the case for incident P waves in Fig.21(b).

For SV waves with under-critical incident angles, the trend of displacement amplitudes for increasing frequency and the trend for those along different directions has nothing special, compared with the case for incident P waves. Thus, analysis will not be presented for this part.

7.3 Time domain responses of half-space for P-wave incidence

The purpose of the following two sections is to analyze the responses of acceleration for oblique incident seismic waves in time domain using the Taft Earthquake mentioned in Subsection 4.1.

Assume that the P-wave vertical acceleration time history given in Fig.7(b) for point OP-1 (see Fig.9) remains unchanged. But three different incident angles, 15°, 50°, 85° will be considered in the analysis. All the other parameters of the model are identical to those used in Section 6.3. The horizontal and vertical ground acceleration time histories \left(a_x,a_y\right), respectively, at different ground observation points were plotted in Fig.26 and Fig.27. Obviously, for both the horizontal and vertical acceleration time histories \left(a_x,a_y\right), the increase in incident angle has resulted in more pronounced phase delay becomes, as well as in amplitude attenuation. Besides, one observes that the horizontal ground acceleration time history \left(a_x\right) at point OP-1 varies for different incident angles, but the vertical ground acceleration time history remains unchanged.

Subsequently, to investigate the influence of incident angle on time domain response, the maximum values of horizontal acceleration time histories \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}}\right) on the free surface and along the vertical central line of the soil model caused by P waves with different incident angles were analyzed for the vertical acceleration time history \left(a_y\right) imposed at point OP-1. As can be observed from Fig.28, the contribution of seismic waves to the horizontal acceleration responses \left(a_x\right) increases with increase in incident angle under the P-waves. In this regard, one can also infer that the vertical acceleration responses \left(a_y\right) decrease with increasing incident angle of incidence P waves.

To visualize the trend of acceleration response along the horizontal and depth directions, the maximum values of the horizontal and vertical acceleration time histories \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{a}\mathrm{n}\mathrm{d}{a}_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) along the x- and y-axes, respectively, under incident P waves for the incident angle {\theta }_{\mathrm{P}\mathrm{i}} = 50° were studied and plotted in Fig.29(a) and 29(b). Besides, the horizontal axis in Fig.29(b) represents the maximum values of acceleration time histories, which will also be adopted in the next sections. From Fig.29, one observes that the maximum values of the horizontal and vertical acceleration time histories \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{a}\mathrm{n}\mathrm{d}{a}_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) uniformly attenuate in the form of fluctuation along the horizontal direction, but irregularly attenuate along the depth direction. Moreover, the maximum horizontal and vertical acceleration responses \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{a}\mathrm{n}\mathrm{d}{a}_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) attenuate more sharply near the ground along the depth direction, but not so serious for greater depth.

7.4 Time domain responses of half-space for SV-wave incidence

Similar to the analysis for the P-wave incidence, the horizontal and vertical ground acceleration time histories \left(a_x,a_y\right), respectively, at the surface observation points caused by SV waves with different incident angles were plotted in Fig.30 and Fig.31, subjected to the horizontal acceleration time history shown in Fig.7(a) for point OP-1 given in Fig.9(a). As can be seen, except for the characteristics of ground responses similar to those caused by incident P waves, the phase delay and amplitude attenuation on the free surface along the x-axis are more obvious, by comparing the peak deviation in time history curves under two types of waves with the same incident angle.

Next, the maximum values of the vertical acceleration time histories \left(a_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) on the free surface and vertical central line of the soil model, respectively, were plotted in Fig.32(a) and Fig.32(b), subjected to the horizontal acceleration time history \left(a_x\right) of point OP-1 given in Fig.9(a). Of interest is that for increasing incident angle, the maximum vertical acceleration response \left(a_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) reaches a peak at the critical angle and then drops to values higher than those prior to the critical angle. The fact that the incidence of SV waves over critical angle can result in generally larger peaks in the vertical responses should be given more attention in practice.

Finally, the maximum values of both the horizontal and vertical acceleration time histories \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{a}\mathrm{n}\mathrm{d}{a}_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) along the horizontal and depth directions, respectively, for SV waves with incident angle {\theta }_{\mathrm{S}\mathrm{i}} = 50° were plotted in Fig.33(a) and Fig.33(b). Clearly, the trends of maximum values of acceleration time history \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{a}\mathrm{n}\mathrm{d}{a}_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) along the horizontal direction are similar to those for incident P waves, as was revealed in Fig.33(a), but with more fluctuations in the responses. All the maximum vertical acceleration responses \left(a_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) are much larger than the horizontal ones \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}}\right), which once more confirms the dramatic amplification effect of the vertical response induced by incidence of SV waves over the critical angle, as shown in Fig.13, Fig.22, Fig.23, Fig.25(a), and Fig.32.

However, the trends of maximum values of acceleration time history \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{a}\mathrm{n}\mathrm{d}{a}_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) in two directions show obvious differences along the depth direction in Fig.33(b). That is, the maximum vertical acceleration response \left(a_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right) reach the peak on the ground, but attenuate obviously with increasing depth, and fluctuate for greater depth, while those of the horizontal ones \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}}\right) tend to increase generally along the depth direction. Moreover, significant phase difference exists between the overall horizontal acceleration response \left(a_{x\mathrm{m}\mathrm{a}\mathrm{x}}\right) and vertical ones \left(a_{y\mathrm{m}\mathrm{a}\mathrm{x}}\right), which is consistent with the trend of the response in frequency domain revealed in Fig.25(b). It can be explained that the absence of reflected P waves and the occurrence of surface waves induce the nonsynchronous motion of particles in the x- and y-axes directions. But for the case where SV waves incident with angles less than the critical one, the trends of maximum acceleration time history are similar to those for incident P waves.

8 Conclusions

This research presents the fundamental theory for the 2D dynamic analysis of a viscoelastic half-space response to oblique incident seismic waves (P or SV waves). Based on the exact solutions in frequency-space domain for the free-field displacements and stresses and the earthquake spectrum, the equivalent seismic forces imposed on the near-field boundary are calculated. Subsequently, using the (asymmetric) finite/infinite element model, the seismic responses of the half-space in frequency and time domains are solved by using the MATLAB software. The reliability of the proposed approach has been verified against Wolf’s and Luco and De Barros’ solutions, as well as for inversely computed ground motions.

Based on the theory, model, material parameters, and earthquake data adopted, some conclusions for the half-space under the 2D oblique incident seismic waves are drawn below.

1) Frequency domain: (a) for both incident P and SV waves, the responses of displacement at larger depth are more sensitive to the change in incident angle; (b) for SV waves with critical incident angle, abrupt change in ground or underground displacement amplitudes occurs.

2) Time domain: (a) the increase in incident angle results in more pronounced phase delay of the acceleration time history, as well as in amplitude attenuation, especially for SV waves; (b) for incident P waves, the maximum horizontal acceleration response of the half-space increases with increasing incident angle, but conversely the vertical response; (c) for incident SV waves, the maximum vertical acceleration response shows a distinct peak value at the critical angle and then drop to values higher than those prior to the critical angle.

3) Spatial irregularity in frequency domain: (a) for both incident P and SV waves, the horizontal and vertical displacements monotonously decrease along the horizontal direction; (b) as for the depth direction, the horizontal and vertical displacements fluctuate synchronously for incident P waves, but there exists half a period of fluctuation for SV waves with over-critical incident angles, both with little attenuation.

4) Spatial irregularity in time domain: (a) for both incident P and SV waves, the maximum accelerations in two directions attenuate and fluctuate along the horizontal direction; (b) for the depth direction, the maximum horizontal and vertical accelerations attenuate sharply near the ground for incident P waves, while significant phase difference exists in the overall acceleration responses for two directions for SV waves with over-critical incident angles.

5) Effect of frequency for both P and SV waves: (a) the seismic waves with low-frequency mainly dominate the ground responses for the earthquake considered; (b) the underground responses fluctuate for increasing frequency with little attenuation owing to the soil’s resonance effect.

Further studies will be conducted to deal with the dynamic responses of the viscoelastic half-space for the 3D oblique incidence of seismic waves, including those of the underground structures.

References

[1] SigakiTKiyohara KSonoYKinositaDMasaoT TamuraRYoshimura CUgataT. Estimation of earthquake motion incident angle at rock site. In: Proceedings of 12th World Conference on Earthquake Engineering. In: Proceedings of 12th World Conference on Earthquake Engineering, 2000, 1–8

[2] EwingW MJardetzky W SPressF. Elastic Waves in Layered Media. New York: McGraw-Hill, 1957

[3] GraffK F. Wave Motion in Elastic Solids. New York: Dover Publications, Inc, 1975

[4] WolfJ P. Dynamic Soil-structure Interaction. Englewood Cliffs: Prentice-Hall, 1985

[5] Lee V W, Karl J. Diffraction of SV waves by underground, circular, cylindrical cavities. Soil Dynamics and Earthquake Engineering, 1992, 11(8): 445–456

[6] Yuan X, Liao Z. Scattering of plane SH waves by a cylindrical alluvial valley of circular-arc cross-section. Earthquake Engineering & Structural Dynamics, 1995, 24(10): 1303–1313

[7] Davis C A, Lee V W, Bardet J P. Transverse response of underground cavities and pipes to incident SV waves. Earthquake Engineering & Structural Dynamics, 2001, 30(3): 383–410

[8] Liang J, Yan L, Lee V W. Scattering of incident plane P waves by a circular-arc canyon with a covering layer. Acta Mechanica Solida Sinica, 2002, 23(4): 397–411

[9] Wong H L. Effect of surface topography on the diffraction of P, SV, and Rayleigh waves. Bulletin of the Seismological Society of America, 1982, 72(4): 1167–1183

[10] Wong K C, Shah A H, Datta S K. Dynamic stresses and displacements in a buried tunnel. Journal of Engineering Mechanics, 1985, 111(2): 218–234

[11] de Barros F C P, Luco J E. Diffraction of obliquely incident waves by a cylindrical cavity embedded in a layered viscoelastic half-space. Soil Dynamics and Earthquake Engineering, 1993, 12(3): 159–171

[12] Luco J E, de Barros F C P. Dynamic displacements and stresses in the vicinity of a cylindrical cavity embedded in a half-space. Earthquake Engineering & Structural Dynamics, 1994, 23(3): 321–340

[13] Luco J E, de Barros F C P. Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space. I: Formulation. Earthquake Engineering & Structural Dynamics, 1994, 23(5): 553–567

[14] de Barros F C P, Luco J E. Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space. II. Validation and numerical results. Earthquake Engineering & Structural Dynamics, 1994, 23(5): 569–580

[15] Luco J E, de Barros F C P. Three-dimensional response of a layered cylindrical valley embedded in a layered half-space. Earthquake Engineering & Structural Dynamics, 1995, 24(1): 109–125

[16] Stamos A A, Beskos D E. Dynamic analysis of large 3-D underground structures by the BEM. Earthquake Engineering & Structural Dynamics, 1995, 24(6): 917–934

[17] Stamos A A, Beskos D E. 3-D seismic response analysis of long lined tunnels in half-space. Soil Dynamics and Earthquake Engineering, 1996, 15(2): 111–118

[18] Mostafa Shaaban A, Rashed Y F. A coupled BEM-stiffness matrix approach for analysis of shear deformable plates on elastic half space. Engineering Analysis with Boundary Elements, 2013, 37(4): 699–707

[19] Ba Z, Yin X. Wave scattering of complex local site in a layered half-space by using a multidomain IBEM: incident plane SH waves. Geophysical Journal International, 2016, 205(3): 1382–1405

[20] Ba Z, Lee V W, Liang J, Yan Y. Scattering of plane qP-and qSV-waves by a canyon in a multi-layered transversely isotropic half-space. Soil Dynamics and Earthquake Engineering, 2017, 98: 120–140

[21] Ba Z N, Zhang E W, Liang J W, Lu Y, Wu M T. Two-dimensional scattering of plane waves by irregularities in a multi-layered transversely isotropic saturated half-space. Engineering Analysis with Boundary Elements, 2020, 118: 169–187

[22] Ba Z, Fu J, Liu Y, Lee V W, Wang Y. Scattering of elastic spherical P, SV, and SH waves by three-dimensional hill in a layered half-space. Soil Dynamics and Earthquake Engineering, 2021, 147: 06545

[23] Liang J, Wang Y, Ba Z, Zhong H. A special indirect boundary element method for seismic response of a 3D canyon in a saturated layered half-space subjected to obliquely incident plane waves. Engineering Analysis with Boundary Elements, 2021, 132: 182–201

[24] Panji M, Mojtabazadeh-Hasanlouei S. On subsurface box-shaped lined tunnel under incident SH-wave propagation. Frontiers of Structural and Civil Engineering, 2021, 15(4): 948–960

[25] Mostafa Shaaban A, Anitescu C, Atroshchenko E, Rabczuk T. Isogeometric boundary element analysis and shape optimization by PSO for 3D axi-symmetric high frequency Helmholtz acoustic problems. Journal of Sound and Vibration, 2020, 486: 115598

[26] Mostafa Shaaban A, Anitescu C, Atroshchenko E, Rabczuk T. Shape optimization by conventional and extended isogeometric boundary element method with PSO for two-dimensional Helmholtz acoustic problems. Engineering Analysis with Boundary Elements, 2020, 113: 156–169

[27] Mostafa Shaaban A, Anitescu C, Atroshchenko E, Rabczuk T. 3D isogeometric boundary element analysis and structural shape optimization for Helmholtz acoustic scattering problems. Computational Methods in Applied Mathematics, 2021, 384: 113950

[28] Mostafa Shaaban A, Anitescu C, Atroshchenko E, Alajlan N, Rabczuk T. Numerical investigations with extended isogeometric boundary element analysis (XIBEM) for direct and inverse Helmholtz acoustic problems. Engineering Analysis with Boundary Elements, 2022, 143: 535–546

[29] Mostafa Shaaban A, Anitescu C, Atroshchenko E, Rabczuk T. An isogeometric Burton-Miller method for the transmission loss optimization with application to mufflers with internal extended tubes. Applied Acoustics, 2022, 185: 108410

[30] Komatitsch D, Vilotte J P, Vai R, Castillo-Covarrubias J M, Sánchez-Sesma F J. The spectral element method for elastic wave equations—application to 2-D and 3-D seismic problems. International Journal for Numerical Methods in Engineering, 1999, 45(9): 1139–1164

[31] Liu J, Du Y, Du X, Wang Z, Wu J. 3D viscous-spring artificial boundary in time domain. Earthquake Engineering and Engineering Vibration, 2006, 5(1): 93–102

[32] Liu J, Gu Y, Li B, Wang Y. An efficient method for the dynamic interaction of open structure-foundation systems. Frontiers of Structural and Civil Engineering, 2007, 1(3): 340–345

[33] Du J, Lin G. Improved numerical method for time domain dynamic structure-foundation interaction analysis based on scaled boundary finite element method. Frontiers of Structural and Civil Engineering, 2008, 2(4): 336–342

[34] Du X, Zhao M. Stability and identification for rational approximation of frequency response function of unbounded soil. Earthquake Engineering & Structural Dynamics, 2010, 39(2): 165–186

[35] Hatzigeorgiou G D, Beskos D E. Soil–structure interaction effects on seismic inelastic analysis of 3-D tunnels. Soil Dynamics and Earthquake Engineering, 2010, 30(9): 851–861

[36] Yu H T, Yuan Y, Bobet A. Multiscale method for long tunnels subjected to seismic loading. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37(4): 374–398

[37] Zhao M, Yin H, Du X, Liu J, Liang L. 1D finite element artificial boundary method for layered half-space site response from obliquely incident earthquake. Earthquakes and Structures, 2015, 9(1): 173–194

[38] Mossaiby F, Shojaei A, Boroomand B, Zaccariotto M, Galvanetto U. Local Dirichlet-type absorbing boundary conditions for transient elastic wave propagation problems. Computational Methods in Applied Mathematics, 2020, 362: 112856

[39] Shojaei A, Hermann A, Seleson P, Cyron C J. Dirichlet absorbing boundary conditions for classical and peridynamic diffusion-type models. Computational Mechanics, 2020, 66(4): 773–793

[40] Huang J Q, Du X, Jin L, Zhao M. Impact of incident angles of P waves on the dynamic responses of long lined tunnels. Earthquake Engineering & Structural Dynamics, 2016, 45(15): 2435–2454

[41] Yan L, Haider A, Li P, Song E. A numerical study on the transverse seismic response of lined circular tunnels under obliquely incident asynchronous P and SV waves. Tunnelling and Underground Space Technology, 2020, 97: 103235

[42] UnglessR F. Infinite finite element. Dissertation for the Master’s Degree. Vancouver: University of British Columbia., 1973

[43] Bettess P. Infinite elements. International Journal for Numerical Methods in Engineering, 1977, 11(1): 53–64

[44] Zhao C, Valliappan S, Wang Y C. A numerical model for wave scattering problems in infinite media due to p-and sv-wave incidences. International Journal for Numerical Methods in Engineering, 1992, 33(8): 1661–1682

[45] Zhao C, Valliappan S. Seismic wave scattering effects under different canyon topographic and geological conditions. Soil Dynamics and Earthquake Engineering, 1993, 12(3): 129–143

[46] Zhao C, Valliappan S. Incident P and SV wave scattering effects under different canyon topographic and geological conditions. International Journal for Numerical and Analytical Methods in Geomechanics, 1993, 17(2): 73–94

[47] Zhao C, Valliappan S. A efficient wave input procedure for infinite media. Communications in Numerical Methods in Engineering, 1993, 9(5): 407–415

[48] Yang Y B, Kuo S R, Hung H H. Frequency-independent infinite element for analyzing semi-infinite problems. International Journal for Numerical Methods in Engineering, 1996, 39(20): 3553–3569

[49] Yun C B, Kim D K, Kim J M. Analytical frequency-dependent infinite elements for soil–structure interaction analysis in two-dimensional medium. Engineering Structures, 2000, 22(3): 258–271

[50] Kim D K, Yun C B. Time-domain soil-structure interaction analysis in two-dimensional medium based on analytical frequency-dependent infinite elements. International Journal for Numerical Methods in Engineering, 2000, 47(7): 1241–1261

[51] Wang G, Chen L, Song C. Finite–infinite element for dynamic analysis of axisymmetrically saturated composite foundations. International Journal for Numerical Methods in Engineering, 2006, 67(7): 916–932

[52] Kouroussis G, Verlinden O, Conti C. Ground propagation of vibrations from railway vehicles using a finite/infinite-element model of the soil. Proceedings of the Institution of Mechanical Engineers. Part F, Journal of Rail and Rapid Transit, 2009, 223(4): 405–413

[53] Yang Y B, Hung H H, Lin K C, Cheng K W. Dynamic response of elastic half-space with cavity subjected to P and SV waves by finite/infinite element approach. International Journal of Structural Stability and Dynamics, 2015, 15(7): 1540009

[54] Lin K C, Hung H H, Yang J P, Yang Y B. Seismic analysis of underground tunnels by the 2. 5D finite/infinite element approach. Soil Dynamics and Earthquake Engineering, 2016, 85: 31–43

[55] Lin S Y, Hung H H, Yang J P, Yang Y B. Seismic analysis of twin tunnels by a finite/infinite element approach. International Journal of Geomechanics, 2017, 17(9): 04017060

[56] Hung H H, Yang Y B. Elastic waves in visco-elastic half-space generated by various vehicle loads. Soil Dynamics and Earthquake Engineering, 2001, 21(1): 1–17

[57] SeedH BIdriss I M. Soil Moduli and Damping Factors for Dynamic Response Analysis. Report No EERC 70-10. 1970

[58] Rayleigh L. On waves propagated along the plane surface of an elastic solid. In: Proceedings of London Mathematical Society, 1885, 1(1): 4–11

[59] Bettess P, Zienkiewicz O C. Diffraction and refraction of surface waves using finite and infinite elements. International Journal for Numerical Methods in Engineering, 1977, 11(8): 1271–1290

RIGHTS & PERMISSIONS

Higher Education Press