A numerical framework for underground structures in layered ground under inclined P-SV waves using stiffness matrix and domain reduction methods

Abstract

A numerical framework was proposed for the seismic analysis of underground structures in layered ground under inclined P-SV waves. The free-field responses are first obtained using the stiffness matrix method based on plane-wave assumptions. Then, the domain reduction method was employed to reproduce the wavefield in the numerical model of the soil–structure system. The proposed numerical framework was verified by providing comparisons with analytical solutions for cases involving free-field responses of homogeneous ground, layered ground, and pressure-dependent heterogeneous ground, as well as for an example of a soil–structure interaction simulation. Compared with the viscous and viscous-spring boundary methods adopted in previous studies, the proposed framework exhibits the advantage of incorporating oblique incident waves in a nonlinear heterogeneous ground. Numerical results show that SV-waves are more destructive to underground structures than P-waves, and the responses of underground structures are significantly affected by the incident angles.

Graphical abstract

Keywords

underground structures / seismic response / stiffness matrix method / domain reduction method / P-SV waves

1 Introduction

The response of underground structures subjected to seismic loading is an important topic in geotechnical earthquake engineering. With the rapid development of computer technology, numerical simulation presents a promising approach for accomplishing seismic analysis of underground structures. Generally, numerical simulations of seismic responses of underground structures should consider critical issues such as the input and radiation of earthquake motions at the truncated boundary of numerical models, wave propagation in the near ground, soil–structure interaction, and dynamic behavior of underground structures [1]. With respect to the latter three issues, considerable research has been conducted in recent years via benchmark tests [2,3] or numerical validation [4] to evaluate the capabilities of the available numerical models. However, the first issue, i.e., the input and radiation of earthquake motions, is beyond the scope of extant studies, in which only vertically incident motions (one-dimensional problem) are considered.

To address this issue, some researchers [1,5] developed comprehensive numerical models to simultaneously consider fault ruptures around the hypocenter, seismic wave propagation in heterogeneous ground, and engineering structures in a particular site as schematically illustrated in Fig.1(a). These advanced models can break the explicit or implicit assumptions including earthquake source and site condition. However, it is difficult to incorporate high-frequency components of earthquake motions and strong nonlinearity of near-surface soils in these models. Moreover, the exceedingly heavy computational burden of such a large-scale region also prohibits the use of this scheme as a routine procedure for analyzing the seismic responses of most underground structures.

A more common and economical practice involves introducing the assumption of plane waves, which is widely and implicitly used in current numerical models for seismic analysis of underground structures because the distance between the source and most underground structures is sufficiently large. Based on Huygens’ principle, the seismic wavefield in the local model, as shown in Fig.1(b), can be recovered if the displacements at the truncated boundary are accurately replicated. The earthquake source and wave propagation history can be partially represented by the characteristics of the input earthquake motion at the truncated boundary of the numerical model. When compared with the former scheme in Fig.1(a), this method exhibits the advantage of significantly lower computational cost as well as the incorporation of nonlinear behaviors of soil and structures. Løkke and Chopra [6] used viscous boundaries [7] to input earthquake motions for a dam-water-foundation interaction system and obtained the wavefield via one-dimensional deconvolution. Zhao et al. [8] used viscous boundaries to analyze tunnel responses under obliquely incident P-SV waves. They used the wave potential method to obtain the free-field response of layered ground. Li and Song [9] employed viscous-spring boundaries and corresponding equivalent nodal force methods to analyze the longitudinal responses of tunnels under inclined P, SV, and SH waves. The free-field analysis was performed based on a soil column with viscous boundaries at the bottom. Huang et al. [10] analyzed the seismic responses of tunnels under P, SV, and SH waves using the same viscous-spring boundary method. An analytical solution for homogeneous ground under inclined waves was used to obtain the free-field responses. Subsequently, the same method was applied by Sun et al. [11] for hydraulic tunnels. The paradigms used in these analyses were similar. First, the free field is obtained, and then the wavefield is recovered in the numerical models via the viscous boundary or viscous-spring boundary. However, viscous and viscous-spring boundaries are local artificial boundaries, which inevitably suffer from long-term instability issues [12] and other limitations of artificial boundaries. Hence, to a certain extent, most of the aforementioned numerical simulations consider only homogeneous ground. Additionally, the two types of boundaries were originally proposed for internal-source problems; therefore, they are inherently inconsistent with the seismic analysis of underground structures, which is a typical external-source problem.

To remedy the aforementioned flaws, the two steps of the paradigm, namely free-field analysis and wavefield reconstruction, should be improved. Different analytical methods are available for the free-field response of layered ground, including the transfer matrix method [13], thin-layer method [14], transmission and reflection method [15], and stiffness matrix method [16]. The stiffness matrix method, which was proposed by Kausel and Roësset [16], yielded symmetric, stable, and well-conditioned matrices with negative exponential terms. The exact stiffness matrices associated with displacements and stresses were formulated for each soil layer, and the global equilibrium equation was then assembled. Hence, this similarity to the finite-element method leads to higher compatibility with numerical simulations. Therefore, the stiffness matrix method was adopted to obtain the wavefield for the inclined P-SV waves. For wavefield reconstruction in numerical models, domain reduction method (DRM) is more advanced than viscous and viscous-spring methods. DRM was proposed by Bielak et al. [17] to analyze local site effects in seismology. It is deduced based on the finite element method for external-source problems, in which the ground near structures can be considered as elastic or elastoplastic, homogeneous, or heterogeneous. Regarding analysis under inclined motions, Zhang et al. [18] used this method to analyze the responses of a rectangular underground structure under obliquely incident SV waves in homogeneous ground. They further developed their method for a layered ground [19]. However, the potential method employed to obtain the free-field responses has intrinsic limitations. This method introduces positive exponential terms, which can easily lead to arithmetic overflow [16]. Compared with the potential method, the stiffness matrix method is stable and well-behaved, and thereby, facilitates the free-field responses of heterogeneous ground as illustrated in Subsection 4.3.

The objective of this study involves proposing a novel numerical framework for the seismic analysis of underground structures in layered ground under obliquely incident P-SV waves. The free-field responses of the layered ground are first obtained using the stiffness matrix method in Section 2, and then DRM is used to recover the wavefield in the numerical model as illustrated in Section 3. The proposed framework outperforms the methods available in the literature because it is stable and applicable to nonlinear heterogeneous ground subjected to inclined plane waves. The verification of the proposed framework is presented in Section 4. Then, examples of free-field responses and soil–structure interactions are provided in Section 5. Finally, the conclusions of the study are summarized in Section 6.

2 Free field responses of layered ground

The prerequisite of the free-field analysis is the stiffness matrix of a homogeneous soil layer, which is presented in Electronic Supplementary Material. Then, based on the continuity and equilibrium conditions between adjacent layers, the global dynamic stiffness matrix of each layered ground, as shown in Fig.2, can be assembled; thus, the global dynamic equilibrium equation can be obtained as follows:

(1)\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ ⋮\\ \tilde{\mathit{L}}\end{array}\right]=\left[\begin{array}{cccccc}{\mathit{K}}_{11}& {\mathit{K}}_{12}& 0& 0& \cdots & 0\\ {\mathit{K}}_{21}& {\mathit{K}}_{22}& {\mathit{K}}_{23}& 0& \cdots & 0\\ 0& {\mathit{K}}_{32}& {\mathit{K}}_{33}& {\mathit{K}}_{34}& \cdots & 0\\ 0& 0& {\mathit{K}}_{43}& {\mathit{K}}_{44}& \ddots & \cdots \\ ⋮& ⋮& ⋮& \ddots & \ddots & {\mathit{K}}_{n-1,n}\\ 0& 0& 0& \cdots & {\mathit{K}}_{n,n-1}& {\mathit{K}}_{nn}+{\mathit{K}}_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}}\end{array}\right]\left[\begin{array}{c}{\tilde{\mathit{U}}}_1\\ {\tilde{\mathit{U}}}_2\\ {\tilde{\mathit{U}}}_3\\ {\tilde{\mathit{U}}}_4\\ ⋮\\ {\tilde{\mathit{U}}}_n\end{array}\right],

where {\mathit{K}}_{ij} denotes a 2\times 2 submatrix of the stiffness matrices of the soil layers. Furthermore, {\tilde{\mathit{U}}}_i denotes the displacement at every soil layer interface in the frequency–wavenumber domain, and \tilde{\mathit{L}} denotes the load due to the inclined P-SV waves at the top surface of the underlying half-space, which is as follows:

(2)\tilde{\mathit{L}}=A_{\mathrm{i}\mathrm{n}\mathrm{p}}\frac{2k\mu (1-s^2)}{1-ps}\left[\begin{array}{cc}p& 0\\ 0& s\end{array}\right]\left(\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ -\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{P}\text{-}\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e},

(3)\tilde{\mathit{L}}=A_{\mathrm{i}\mathrm{n}\mathrm{p}}\frac{2k\mu (1-s^2)}{1-ps}\left[\begin{array}{cc}p& 0\\ 0& s\end{array}\right]\left(\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ \mathrm{i}\mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{V}\text{-}\mathrm{w}\mathrm{a}\mathrm{v}\mathrm{e},

where A_{\mathrm{i}\mathrm{n}\mathrm{p}} denotes the Fourier amplitude of the input displacement. Viscous damping can be incorporated by the well-known correspondence principle, namely by replacing the Lamé constants with complex constants.

(4)\begin{array}{rl}& {\lambda }^{\ast }=\lambda \left(1+2\xi \mathrm{i}\right)\\ & {\mu }^{\ast }=\mu (1+2\xi \mathrm{i}),\end{array}

where \xi denotes the damping ratio. The free-field responses in the time–spatial domain can be obtained using the inverse double Fourier transform used in the DRM step. It should be noted that complex damping is used in free-field analysis, which differs from the widely used Rayleigh damping in numerical models.

3 Wavefield reconstruction in numerical models

The DRM is adopted to reconstruct the wavefield in the numerical model for the seismic analysis of underground structures. In DRM [17], as shown in Fig.3, the outer boundary {\hat{\mathrm{\Gamma }}}^{+} truncates the numerical model from an infinite domain. The entire numerical model is divided by boundary Γ into two parts: {\mathrm{\Omega }}^{-} and {\mathrm{\Omega }}^{+}. The interior domain {\mathrm{\Omega }}^{-} contains the underground structure of interest, capturing its seismic responses as well as the soil–structure interaction. The exterior domain {\mathrm{\Omega }}^{+} is an auxiliary region that applies seismic excitations and dissipates possible outgoing waves from the interior domain. An auxiliary boundary exists for applying seismic excitations in the exterior domain, {\mathrm{\Omega }}^{+}. The equivalent nodal forces P^{\text{eff}} are imposed on the nodes in a single-element layer enclosed by the boundaries Γ and {\mathrm{\Gamma }}_{\mathrm{e}} as follows:

(5)\left[\begin{array}{c}{\mathit{P}}_{\mathrm{b}}^{\mathrm{e}\mathrm{f}\mathrm{f}}\\ {\mathit{P}}_{\mathrm{e}}^{\mathrm{e}\mathrm{f}\mathrm{f}}\end{array}\right]=\left[\begin{array}{c}-{\mathit{M}}_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}{{\ddot{\mathit{u}}}_{\mathrm{e}}}^{\mathrm{f}}-{\mathit{K}}_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}{\mathit{u}}_{\mathrm{e}}^{\mathrm{f}}\\ {\mathit{M}}_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}{{\ddot{\mathit{u}}}_{\mathrm{b}}}^{\mathrm{f}}+{\mathit{K}}_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}{\mathit{u}}_{\mathrm{b}}^{\mathrm{f}}\end{array}\right].

The detailed process for deducing the equivalent nodal force is presented in Electronic Supplementary Material. Given that the possible outgoing waves from the interior domain due to the underground structures are expected to be dissipated by the region enclosed by {\mathrm{\Gamma }}_{\mathrm{e}} and {\hat{\mathrm{\Gamma }}}^{+}, additional Rayleigh damping reduces the computational cost [20]. The dimensions of the region and corresponding Rayleigh damping are related to the plasticity of the interior domain, underground structure, seismic motion. Lower damping ratios can be used when the region is large and vice versa. The dimensions of the region and damping imposed can be selected through preliminary simulations in specific applications. The method in this study was implemented in Abaqus [21] and OpenSees [22] using Python script. In the implementation, the thicknesses of the soil layers in the stiffness matrix method are consistent with the element dimensions of the numerical model along the depth to incorporate the variation of soil parameters in the vertical direction.

4 Verification

4.1 Homogeneous ground

Consider a homogeneous elastic ground with a Young’s modulus E_s = 250 MPa, Poisson’s ratio v_s = 0.25, and density {\rho }_{\mathrm{e}} = 2000 kg/m³. A Ricker wavelet with a central frequency of 3 Hz and an amplitude of 1 m, as shown in Fig.4, is input from a depth of 200 m as P-waves. The incidence angle was 75°. The width of the interior domain was 200 m. The observation points were on the medial axis of the numerical model (x = 0). The seismogram comparison between the analytical and numerical solutions is shown in Fig.5, in which the results are magnified five times for better visualization. This consistency indicated that the analytical wavefield is fully reconstructed in the numerical model. The wave propagation trace is clear: the motion is inputted upward and reflected downward at the ground surface. The displacement responses at points (0,–40) and (0,–140) are compared in Fig.6. The time histories and spectra of the analytical and numerical solutions coincided, and the accuracy of the proposed method was verified. The displacement contours of the numerical model at different moments are shown in Fig.7, where the obliquely incident plane wavefront is shown. The reflection occurs at 1.4–1.8 s, and then the reflected wave propagates downward. It can be observed that the plane wavefronts of the reflected wave are not as clear as those of the incident wave. This is due to the fact that the incident P-wave simultaneously induces reflected P- and SV-waves.

4.2 Layered ground

Consider a two-layer ground, and each layer is 100 m in thickness. The Young’s moduli of the upper and lower layers were E_s1 = 250 MPa and E_s2 = 500 MPa, respectively. The Poisson’s ratio and density of both layers were v_s = 0.25 and {\rho }_{\mathrm{s}}= 2000 kg/m³, respectively. The Ricker wavelet that is input from a depth of 200 m is similar to that of the SV-waves. The incidence angle was 80°. Analytical and numerical seismograms are presented in Fig.8. The consistency between the analytical and numerical solutions again verified that the wavefield was completely reconstructed in the numerical model. The displacements at two positions (0,–60) and (0,–120) are compared in Fig.9. The numerical solutions are consistent with the analytical solutions, which again verify the accuracy of the proposed method. As shown in Fig.9, the displacement responses are more complex than those in the homogeneous ground (Fig.6) owing to reflections and refractions at the interface of the two soil layers.

4.3 Pressure-dependent heterogeneous ground

It is widely known that the maximum shear modulus of granular soil is dependent on void ratios, confining pressures, and other parameters [23,24]. Assuming that the ground is uniform, its maximum shear modulus increases gradually with respect to depth owing to the increasing confining pressure values. By revisiting the homogeneous ground in Subsection 4.1, the following relationship can be considered.

(6)G_0=G_{\mathrm{r}\mathrm{e}\mathrm{f}}{\left(\frac{p}{p_{\mathrm{r}\mathrm{e}\mathrm{f}}}\right)}^{0.5},

where G_{\mathrm{r}\mathrm{e}\mathrm{f}}= 250 MPa and p_{\mathrm{r}\mathrm{e}\mathrm{f}} = 101 kPa corresponds to the reference shear modulus and reference confining pressure, respectively, G_0 and p denote the maximum shear modulus and confining pressure at different depths, respectively, which are controlled by the self-weight of the ground. The same Ricker wavelet was input from a depth of 200 m as SV-waves with an incident angle of 70°. The seismograms obtained from the analytical and numerical solutions are compared in Fig.10. The coincidence between them indicates that the method in this study can incorporate the pressure-dependent characteristics of the ground. The displacement contours at different times are shown in Fig.11. It is evident that the numerical model accurately replicates the incidence and reflection of wave propagation. Given the pressure-dependent feature, the shear modulus, as well as the shear wave velocity, of the ground declines from the bottom to the top. Therefore, the wave propagation direction deviates from the original incident direction and gradually approaches the vertical direction, as marked in Fig.11 (1.2 and 1.3 s), because of the well-known Snell’s law.

4.4 Soil–tunnel model

All the aforementioned cases involve free-field analyses, which verify the ability of the proposed method to recover the analytical plane wavefield in numerical models. The soil–tunnel interaction problem was considered in this study. A tunnel with a radius of 3 m and thickness of 0.35 m is considered, which is simulated by beam elements. The density of the tunnel was {\rho }_{\mathrm{t}} = 2500 kg/m³. The Young’s modulus and Poisson’s ratio were E_t = 35.5 GPa and v_t = 0.2, respectively. The tunnel was bonded to the surrounding ground and had a buried depth of 10 m. The ground parameters used in Subsection 4.1 are adopted. The same Ricker wavelet was input vertically as the SV-waves with an incident angle of 90°. To verify the proposed method, an extended model with a depth of 1000 m is built, and the Ricker wavelet is input from the bottom boundary as shown in Fig.12. A repeatable boundary was used at the lateral boundaries. According to the wave propagation theory, the motion reached the top surface after 4.47 s and then reflected downward. Therefore, the simulation time was 8 s to avoid superfluous reflections at the bottom boundary. Thus, the dynamic responses of the tunnel resting on an elastic half-space under vertically incident SV waves were captured. Comparatively, the motion was input from a depth of 100 m using the method proposed in this study. The maximum displacements of the tunnel along the circumferential direction are compared in Fig.13, which verifies the accuracy of the proposed method as well as its applicability to soil–structure interaction problems.

5 Examples

5.1 Free field responses

As introduced in Section 1, viscous and viscous-spring boundaries are widely used in the literature to simulate vertical motions, as well as inclined motions, from underlying elastic bases. The two methods are compared with the method proposed in this study, and their differences are clarified. Considering the one-dimensional problem, SV-waves are inputted vertically from a depth of 200 m as shown in Fig.14. All vertical degrees of freedom were fixed because vertical SV waves were only associated with horizontal displacements. The equivalent nodal forces related to viscous boundaries, which are input at the bottom nodes, can be calculated as follows [7]:

(7)F_{\mathrm{v}}=2Cv\left(t\right),

where C={\rho }_{\mathrm{s}}{V}_{\mathrm{s}}A denotes the dashpot coefficient, {\rho }_{\mathrm{s}} and V_{\mathrm{s}} denote the density and shear wave velocity of the ground, respectively, and A denotes the cross-sectional area of the bottom boundary. Furthermore, v(t) denotes the velocity–time history at the bottom boundary. Correspondingly, the equivalent nodal forces associated with the viscous-spring boundaries are [9] as follows:

(8)F_{\mathrm{v}\mathrm{s}}=Ku\left(t\right)+2Cv\left(t\right),

where K=0.66A{G}_{\mathrm{s}}/R is the spring coefficient, and G_s denotes the shear modulus of the ground. Furthermore, R denotes the distance between the scattering source and viscous-spring boundary. Additionally, u\left(t\right) denotes the displacement–time history at the bottom boundary.

First, the homogeneous ground in Subsection 4.1 is used, and the same Ricker wavelet is input vertically. In this scenario, parameter R = 200 m, in Eq. (8), is designated. The horizontal displacements at a depth of 190 m and at the top surface are compared in Fig.15. Clearly, all the methods obtained satisfactory results and were in good agreement with the analytical solution. Minor discrepancies were observed when the viscous-spring boundary was used.

Then, the two-layer ground in Subsection 4.2 is used, and the motion is still input from a depth of 200 m. In this scenario, parameter R = 100 m in Eq. (8) is considered because there is a reflection at the interface between the two layers. Analogously, the horizontal displacements at a depth of 190 m and top surface are compared in Fig.16. It is clear that the results from the different methods are consistent with the analytical solution. Minor discrepancies occurred at the viscous spring boundary. The aforementioned comparison indicates that the proposed method is in good agreement with the methods available in the literature.

The differences between the proposed method and viscous or viscous spring boundary are clarified here. In the DRM, the single-element thick layer, enclosed by boundaries Γ and {\mathrm{\Gamma }}_{\mathrm{e}}, plays the role of reconstructing the free wavefield based on the applied nodal forces. Simultaneously, the region enclosed by the boundaries {\mathrm{\Gamma }}_{\mathrm{e}} and {\hat{\mathrm{\Gamma }}}^{+} dissipates the possible outgoing waves (the term {\mathit{w}}_{\mathrm{e}} in Eq. (B3)). Correspondingly, the viscous or viscous-spring boundary plays a dual role in reconstructing the free wavefield and dissipating the outgoing waves. Hence, their performance deteriorates when soil nonlinearity is considered in the interior domain. DRM is deduced based on the assumption that the source of excitation is outside the region of interest. The interior region can be elastic, viscoelastic, or elasto-plastic. However, the viscous or viscous-spring boundary method is based on the excitation source inside the near field. Hence, parameter R is included in Eq. (8). Therefore, DRM is more appropriate for the seismic analysis of engineering structures without any additional parameters, and the framework in this study has significant advantages over the viscous and viscous-spring boundary methods. In summary, the numerical framework in this study can facilitate the development of current numerical models used for one-dimensional problems to incorporate wavefields under inclined P-SV waves without further requirements on the ground. The next example demonstrates the framework for the seismic analysis of underground structures with nonlinear soil.

5.2 Soil–structure interaction

A typical rectangular underground structure in a homogeneous and pressure-dependent ground is considered as shown in Fig.17. The burial depth was 10 m. The width and height of the structure were 16 and 6 m, respectively. The thicknesses of the ceiling, floor, and side walls were 1 m while that of the central column was 0.4 m. The structure was simulated by elastic beam elements in OpenSees [22] with a Young’s modulus of 35.5 GPa. The ground was simulated by plane strain elements and the PDMY02 model [25]. The parameters suggested by the developer are listed in Tab.1. For simplicity, the interface between the soil and underground structure is assumed as fully bonded. It should be noted that other complex soil–structure interaction models can be used without any challenges. The Ricker wavelet in Fig.4 is still used in this interaction model and is input from a depth of 50 m. Its displacement amplitude was scaled to 0.002 m, corresponding to an acceleration amplitude of 0.1 g. The Newmark method [26] was used to integrate the dynamic response with parameters \gamma = 0.6 and \beta = 0.3025, which introduce small numerical damping for better convergence.

First, we input the Ricker wavelet as P- and SV-waves with the same incident angle of 75° to compare their differences. The deformation of the structure is shown in Fig.18. The horizontal drift ratio is defined as the ratio of the horizontal displacement difference between the ceiling and floor to the height. The vertical drift ratio is defined as the ratio of the vertical displacement difference of the side walls to the width, which is seldom considered in previous studies. It can be observed that SV-waves induce larger deformations, in terms of amplitudes and residues, than P-waves, whereas P-waves trigger the responses of the structure earlier owing to higher wave velocities. Additionally, it can be observed that the structure rotates toward the right and floor direction first, and the rotation direction is reversed when the phase of motion changes. This indicates that the seismic responses of underground structures under inclined motion differ significantly from those under one-dimensional vertical shear waves. The internal forces on the central column at the top are shown in Fig.19. Furthermore, P-waves excite larger axial forces than SV waves, whereas SV waves excite larger shear forces and bending moments. Additionally, SV waves induce larger permanent internal forces than P-waves. This indicates that SV-waves can excite higher responses of the structure than P-waves, and thus, SV-waves are more destructive.

Then, we input the Ricker wavelet as SV waves with different incident angles of 60°, 70°, 80°, and 90° to investigate the effect of incident angles. Analogously, the deformation of the structure and internal forces of the central column at the top are plotted in Fig.20 and Fig.21, respectively. This clearly suggests that the responses of the structure are affected by the incident angles. The structural deformation amplitudes increased with the incident angle. Larger incident angles excited the responses of the structure more quickly. The shear force and bending moment of the central column exhibited the same tendency as the deformation of the structure. However, the axial force amplitudes decreased and the permanent values increased as incident angles increased. This implies that the assumption of one-dimensional problems, namely the incident angle of 90°, can lead to the underestimating of the axial force of the central column. It is important to note that the numerical framework proposed in this study is compatible with the available elements and constitutive models used in common one-dimensional numerical models. By incorporating inclined waves, the proposed method overcomes the limitations of one-dimensional problem assumptions and provides more accurate numerical results on the seismic responses of underground structures.

6 Conclusions

In this study, a numerical framework for the seismic analysis of underground structures in layered ground under inclined plane waves was proposed and verified. For elastic ground, the method in this study was in good agreement with the viscous and viscous-spring boundary methods. In nonlinear heterogeneous ground, the method outperformed the viscous and viscous-spring boundary methods owing to the advantages of the DRM and stiffness matrix method. The proposed method was used to analyze the seismic responses of rectangular structures. The results showed that SV-waves are more destructive than P-waves, and that the responses of the structure are significantly affected by the incident angles. The deformations of the structure, as well as the shear force and bending moment of the central column, increase with increasing incident angles, whereas the axial force of the central column decreases with increasing incident angles. The proposed method can overcome the limitation of vertical incidence, and thereby, facilitating the seismic analysis of underground structures in nonlinear heterogeneous ground under inclined P-SV waves. Furthermore, it can be used for seismic analysis of other engineering structures. A similar framework for the seismic analysis of underground structures in saturated ground under inclined plane waves will be presented in another follow-up study.

References

[1] McCallen D, Petersson A, Rodgers A, Pitarka A, Miah M, Petrone F, Sjogreen B, Abrahamson N, Tang H. EQSIM—A multidisciplinary framework for fault-to-structure earthquake simulations on exascale computers part I: Computational models and workflow. Earthquake Spectra, 2021, 37(2): 707–735

[2] Bilotta E, Lanzano G, Madabhushi S P G, Silvestri F. A numerical Round Robin on tunnels under seismic actions. Acta Geotechnica, 2014, 9(4): 563–579

[3] Yuan Y, Yang Y S, Zhang S H, Yu H T, Sun J. A benchmark 1 g shaking table test of shallow segmental mini-tunnel in sand. Bulletin of Earthquake Engineering, 2020, 18(11): 5383–5412

[4] Régnier J, Bonilla L F, Bard P Y, Bertrand E, Hollender F, Kawase H, Sicilia D, Arduino P, Amorosi A, Asimaki D, Boldini D, Chen L, Chiaradonna A, DeMartin F, Ebrille M, Elgamal A, Falcone G, Foerster E, Foti S, Garini E, Gazetas G, Gélis C, Ghofrani A, Giannakou A, Gingery J R, Glinsky N, Harmon J, Hashash Y, Iai S, Jeremić B, Kramer S, Kontoe S, Kristek J, Lanzo G, Lernia A, Lopez-Caballero F, Marot M, McAllister G, Diego Mercerat E, Moczo P, Montoya-Noguera S, Musgrove M, Nieto-Ferro A, Pagliaroli A, Pisanò F, Richterova A, Sajana S, Santisi d’Avila M P, Shi J, Silvestri F, Taiebat M, Tropeano G, Verrucci L, Watanabe K. International benchmark on numerical simulations for 1D, nonlinear site response (PRENOLIN): Verification phase based on canonical cases. Bulletin of the Seismological Society of America, 2016, 106(5): 2112–2135

[5] Abell J A, Orbović N, McCallen D B, Jeremic B. Earthquake soil-structure interaction of nuclear power plants, differences in response to 3-D, 3×1-D, and 1-D excitations. Earthquake Engineering & Structural Dynamics, 2018, 47(6): 1478–1495

[6] Løkke A, Chopra A K. Direct finite element method for nonlinear earthquake analysis of concrete dams: Simplification, modeling, and practical application. Earthquake Engineering & Structural Dynamics, 2019, 48(7): 818–842

[7] Lysmer J, Kuhlemeyer R L. Finite dynamic model for infinite media. Journal of the Engineering Mechanics Division, 1969, 95(4): 859–877

[8] Zhao W S, Chen W Z, Yang D S, Tan X J, Gao H, Li C. Earthquake input mechanism for time-domain analysis of tunnels in layered ground subjected to obliquely incident P- and SV-waves. Engineering Structures, 2019, 181: 374–386

[9] Li P, Song E X. Three-dimensional numerical analysis for the longitudinal seismic response of tunnels under an asynchronous wave input. Computers and Geotechnics, 2015, 63: 229–243

[10] Huang J Q, Zhao M, Du X. Non-linear seismic responses of tunnels within normal fault ground under obliquely incident P waves. Tunnelling and Underground Space Technology, 2017, 61: 26–39

[11] Sun B B, Zhang S R, Cui W, Deng M J, Wang C. Nonlinear dynamic response and damage analysis of hydraulic arched tunnels subjected to P waves with arbitrary incoming angles. Computers and Geotechnics, 2020, 118: 103358

[12] Baffet D, Bielak J, Givoli D, Hagstrom T, Rabinovich D. Long-time stable high-order absorbing boundary conditions for elastodynamics. Computer Methods in Applied Mechanics and Engineering, 2012, 241-244: 20–37

[13] Haskell N A. The dispersion of surface waves on multilayered media. Bulletin of the Seismological Society of America, 1953, 43(1): 17–34

[14] Thomson W T. Transmission of elastic waves through a stratified soil medium. Journal of Applied Physics, 1950, 21(2): 89–93

[15] Kausel E. Thin-layer method: Formulation in the time domain. International Journal for Numerical Methods in Engineering, 1994, 37(6): 927–941

[16] Kausel E, Roësset J M. Stiffness matrices for layered soils. Bulletin of the Seismological Society of America, 1981, 71(6): 1743–1761

[17] Bielak J, Loukakis K, Hisada Y, Yoshimura C. Domain reduction method for three-dimensional earthquake modeling in localized regions, part I: Theory. Bulletin of the Seismological Society of America, 2003, 93(2): 817–824

[18] Zhang W, Seylabi E E, Taciroglu E. An ABAQUS toolbox for soil-structure interaction analysis. Computers and Geotechnics, 2019, 114: 103143

[19] Zhang W, Taciroglu E. 3D time-domain nonlinear analysis of soil-structure systems subjected to obliquely incident SV waves in layered soil media. Earthquake Engineering & Structural Dynamics, 2021, 50(8): 2156–2173

[20] Semblat J F, Lenti L, Gandomzadeh A. A simple multi-directional absorbing layer method to simulate elastic wave propagation in unbounded domains. International Journal for Numerical Methods in Engineering, 2011, 85(12): 1543–1563

[21] Abaqus. Version 6.11. Paris: Dassault Systemes Simulia Corporation. 2011

[22] MazzoniSMcKennaFScottM HFenvesG L. OpenSees Command Language Manual. 2006

[23] Miao Y, He H J, Liu H B, Wang S Y. Reproducing ground response using in-situ soil dynamic parameters. Earthquake Engineering & Structural Dynamics, 2022, 51(10): 2449–2465

[24] Wang S Y, Zhuang H Y, Zhang H, He H J, Jiang W P, Yao E L, Ruan B, Wu Y X, Miao Y. Near-surface softening and healing in eastern Honshu associated with the 2011 magnitude-9 Tohoku-Oki Earthquake. Nature Communications, 2021, 12(1): 1215

[25] Elgamal A, Yang Z H, Parra E, Ragheb A. Modeling of cyclic mobility in saturated cohesionless soils. International Journal of Plasticity, 2003, 19(6): 883–905

[26] Newmark N M. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 1959, 85(3): 67–94

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn