An end-to-end 3D seismic simulation of underground structures due to point dislocation source by using an FK-FEM hybrid approach

Abstract

Based on the domain reduction idea and artificial boundary substructure method, this paper proposes an FK-FEM hybrid approach by integrating the advantages of FK and FEM (i.e., FK can efficiently generate high-frequency three translational motion, while FEM has rich elements types and constitutive models). An advantage of this approach is that it realizes the entire process simulation from point dislocation source to underground structure. Compared with the plane wave field input method, the FK-FEM hybrid approach can reflect the spatial variability of seismic motion and the influence of source and propagation path. This approach can provide an effective solution for seismic analysis of underground structures under scenario of earthquake in regions where strong earthquakes may occur but are not recorded, especially when active faults, crustal, and soil parameters are available. Taking Daikai subway station as an example, the seismic response of the underground structure is simulated after verifying the correctness of the approach and the effects of crustal velocity structure and source parameters on the seismic response of Daikai station are discussed. In this example, the influence of velocity structure on the maximum interlayer displacement angle of underground structure is 96.5% and the change of source parameters can lead to the change of structural failure direction.

Graphical abstract

Keywords

source-to-structure simulation / FK-FEM hybrid approach / underground structures / point dislocation source

1 Introduction

Rapid and unreasonable urbanization has brought “city diseases”, such as loss of valuable farm land, environmental pollution, and traffic congestion [1,2]. Underground space development provides an effective solution, such as in the form of subways, road tunnels, utility tunnels, and shopping complexes, which greatly free up the land at surface [3,4]. Due to the constraints of surrounding soil or rock, the seismic performance of underground structures was empirically believed to be superior to that of aboveground structures, thus their seismic designs were rarely accounted [5]. However, in the 1995 Kobe earthquake in Japan, the Daikai subway station and some interval tunnels were seriously damaged, which sounded an alarm for the necessity of seismic design of underground structures. Soon after, scholars worldwide started an upsurge on seismic studies of underground structures [6–10].

The primary methods for seismic analysis of underground structures include the quasi-static analysis method [11–18], dynamic time history analysis method [19,20], and lab test [21,22]. The quasi-static analysis method is straightforward, but ignores the wave effect during seismic wave propagation, such as scattering and reflection. Lab test is an effective, yet costly approach to investigate earthquake damage of underground structures. The dynamic time history analysis method has many advantages, such as clear physical meaning, rich material constitutive models, and the ability to simulate complex models, so it is widely used by researchers [23–25].

It is known that the seismic response of structures not only depends on their own dynamic characteristics, but also on the characteristics of the input seismic motion. For example, different input seismic waves may lead to large deviations in the seismic response of the same structure [26,27]. Therefore, proper input seismic wave is required for seismic evaluations of underground structures [24,28]. In fact, source parameters, wave propagation path, and local site control the characteristics of the input seismic wave field [28–30]. Different site conditions and source parameters will produce distinct seismic wave fields. However, due to the lack of available strong earthquake records in many cities, the prevailing seismic response analysis of underground structures will typically utilize vibration or wave methods to apply seismic load after adjusting the amplitude and spectrum characteristics of existing or artificial seismic waves [31]. In addition, the waveform of the input seismic wave is generally assumed to be plane wave. However, these two methods can not reflect the influence of source parameters and traveling wave effect caused by propagation path and local topography, thus the spatial variability of input seismic cannot be considered. In addition, since the grid size of the structure is much smaller than that of crust, it is difficult to simulate the source-to-underground structures simulation in one model because of computational intractability, which hinders the seismic evaluations of underground structures with source and propagation path included.

The domain reduction idea proposed by Bielak in 2003 solves the problem from source to local topography, which provides a powerful reference for source-to-structures simulations [32–34]. The entire physics-based model is divided into a global domain (containing source and propagation path) and local domain of interest (containing local topography). Firstly, the global domain model is simulated to generate the seismic wave field, which could account for the effect of source and propagation path. Then, the obtained seismic wave field is applied to the local domain model as the input seismic wave. Finally, the source-to-local topography simulation is realized.

Referring to the domain reduction idea, the simulation model is divided into a global domain model (containing source and propagation path) and a local domain model (containing local topography and underground structure). In terms of seismic wave field simulation, the ground motion synthesis method based on Green’s function in frequency wavenumber domain (frequency-wavenumber domain method, i.e., FK method) can produce three translational movements with a broad frequency band from 0.1 to 20 Hz, which has high efficiency and covers the sensitive frequencies of engineering structures [35,36]. Therefore, this paper uses the FK method to simulate the seismic wave propagation in the global domain model with earthquake source included. In terms of structural analysis, the finite element method (FEM) has rich elements types and constitutive models [37,38], and many commercial software (such as ABAQUS, ANSYS, etc.) are available and widely accepted. Therefore, FEM is adopted to simulate the local domain model. In addition, the artificial boundary substructure method [39,40] is used to obtain the equivalent seismic load, which can be applied to the FEM local domain model. Then, the FK-FEM hybrid approach is established to perform the source-to-underground structure simulation, which creates a direct link between the 3D physics-based geophysical simulations and the seismic response analysis of underground structures.

The remainder of this paper is organized as follows. Section 2 introduces the principle and main procedures of the FK-FEM hybrid approach. Section 3 verifies the correctness of the approach by comparing the simulation results from the FK method and FK-FEM hybrid approach. Section 4 verifies the necessity of spatial variability of seismic wave input in underground structure analysis, and discusses in detail the influence of crustal velocity structure and source parameters on the seismic response of underground structure. Section 5 summarizes the main conclusions.

2 Hybrid approach of frequency- wavenumber domain method and finite element method

2.1 Frequency-wavenumber domain method

Displacement response under the action of a dislocation point source can be expressed as Eq. (1) [41],

(1)u_i\left(\mathit{x},t\right)=G\left(\mathit{x},t\right)\ast \mathit{S}\left(t\right){\mathit{M}}_{pq},

where u_i denotes the displacement of the point x along direction i; t denotes time; ‘\ast’ represents the convolution operator; G(x,t) is the dynamical Green’s function under the action of the dislocation point source; S(t) is the source time function; M_pq denotes the pq-th component in the moment tensor, which essentially represents two parallel and identical magnitude forces in p-th and opposite p-th directions distinguished by a perpendicular distance in the q-th direction (p, q = x, y, z).

The main procedures of calculating Green’s function using the FK method is as follows. Firstly, using Fourier Hankel transform, the wave equation in the cylindrical coordinate system of the spatial-temporal domain is transformed into a frequency wavenumber domain. Then, after introducing stress and displacement boundary conditions, the response induced by the source is obtained using the stiffness matrix method. Finally, the response in spatial-temporal domain is obtained by Fourier Hankel inverse transform.

Due to space limitations, derivations of the FK method are not described in this paper. For details on the FK method, please refer to Ref. [42].

2.2 Theory of the hybrid approach

Fig.1 shows that the simulation model is divided into a global domain model (containing source and propagation path) and local domain model (containing local topography and underground structure), according to the domain reduction idea. The substructure of artificial boundary is composed of artificial boundary nodes and internal nodes next to the artificial boundaries.

For the local domain model, the dynamic equilibrium equation could be assembled into block matrix form, as shown in Eq. (2). The dynamic equilibrium equation of the free-field model corresponding to local domain model (i.e., the local domain model that does not contain any structure) is shown in Eq. (3),

(2)\begin{array}{rl}& \left[\begin{array}{ccc}{\mathit{M}}_{\mathrm{I}\mathrm{I}}& {\mathit{M}}_{\mathrm{I}\mathrm{C}}& 0\\ {\mathit{M}}_{\mathrm{C}\mathrm{l}}& {\mathit{M}}_{\mathrm{C}\mathrm{C}}& {\mathit{M}}_{\mathrm{C}\mathrm{B}}\\ 0& {\mathit{M}}_{\mathrm{B}\mathrm{C}}& {\mathit{M}}_{\mathrm{B}\mathrm{B}}\end{array}\right]\left\{\begin{array}{l}{\ddot{\mathit{u}}}_{\mathrm{I}}\\ {\ddot{\mathit{u}}}_{\mathrm{C}}\\ {\ddot{\mathit{u}}}_{\mathrm{B}}\end{array}\right\}+\left[\begin{array}{ccc}{\mathit{C}}_{\mathrm{I}\mathrm{I}}& {\mathit{C}}_{\mathrm{I}\mathrm{C}}& 0\\ {\mathit{C}}_{\mathrm{C}\mathrm{I}}& {\mathit{C}}_{\mathrm{C}\mathrm{C}}& {\mathit{C}}_{\mathrm{C}\mathrm{B}}\\ 0& {\mathit{C}}_{\mathrm{B}\mathrm{C}}& {\mathit{C}}_{\mathrm{B}\mathrm{B}}\end{array}\right]\left\{\begin{array}{l}{\dot{\mathit{u}}}_{\text{I}}\\ {\dot{\mathit{u}}}_{\mathrm{C}}\\ {\dot{\mathit{u}}}_{\mathrm{B}}\end{array}\right\}\\ +& \left[\begin{array}{ccc}{\mathit{K}}_{\mathrm{I}\mathrm{I}}& {\mathit{K}}_{\mathrm{I}\mathrm{C}}& 0\\ {\mathit{K}}_{\mathrm{C}\mathrm{l}}& {\mathit{K}}_{\mathrm{C}\mathrm{C}}& {\mathit{K}}_{\mathrm{C}\mathrm{B}}\\ 0& {\mathit{K}}_{\mathrm{B}\mathrm{C}}& {\mathit{K}}_{\mathrm{B}\mathrm{B}}\end{array}\right]\left\{\begin{array}{l}{\mathit{u}}_{\mathrm{I}}\\ {\mathit{u}}_{\mathrm{C}}\\ {\mathit{u}}_{\mathrm{B}}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ {\mathit{F}}_{\mathrm{B}}\end{array}\right\},\end{array}

(3)\begin{array}{rl}& \left[\begin{array}{ccc}{\mathit{M}}_{\text{II}}^0& {\mathit{M}}_{\text{IC}}^0& 0\\ {\mathit{M}}_{\text{CI}}^0& {\mathit{M}}_{\text{CC}}^0& {\mathit{M}}_{\text{CB}}^0\\ 0& {\mathit{M}}_{\text{BC}}^0& {\mathit{M}}_{\text{BB}}^0\end{array}\right]\left\{\begin{array}{c}{\ddot{\mathit{u}}}_{\text{I}}^0\\ {\ddot{\mathit{u}}}_{\text{C}}^0\\ {\ddot{\mathit{u}}}_{\text{B}}^0\end{array}\right\}+\left[\begin{array}{ccc}{\mathit{C}}_{\text{II}}^0& {\mathit{C}}_{\text{IC}}^0& 0\\ {\mathit{C}}_{\text{CI}}^0& {\mathit{C}}_{\text{CC}}^0& {\mathit{C}}_{\text{CB}}^0\\ 0& {\mathit{C}}_{\text{BC}}^0& {\mathit{C}}_{\text{BB}}^0\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{u}}}_{\text{I}}^0\\ {\dot{\mathit{u}}}_{\text{C}}^0\\ {\dot{\mathit{u}}}_{\text{B}}^0\end{array}\right\}\\ +& \left[\begin{array}{ccc}{\mathit{K}}_{\text{II}}^0& {\mathit{K}}_{\text{IC}}^0& 0\\ {\mathit{K}}_{\text{CI}}^0& {\mathit{K}}_{\text{CC}}^0& {\mathit{K}}_{\text{CB}}^0\\ 0& {\mathit{K}}_{\text{BC}}^0& {\mathit{K}}_{\text{BB}}^0\end{array}\right]\left\{\begin{array}{c}{\mathit{u}}_{\text{I}}^0\\ {\mathit{u}}_{\text{C}}^0\\ {\mathit{u}}_{\text{B}}^0\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ {\mathit{F}}_{\text{B}}\end{array}\right\},\end{array}

where M, C, and K are the mass matrix, damping matrix, and stiffness matrix of the soil-structure model, respectively; F_B is the equivalent input seismic loads imposed on artificial boundary nodes; superscript 0 denotes the free-field model; subscripts B, C, and I denote the artificial boundary nodes, the internal nodes next to the artificial boundary, and remaining internal nodes, respectively.

When the mesh generation of the free-field model at substructure of artificial boundary is consistent with that of the local domain model, Eq. (4) can be obtained,

(4)\begin{array}{lll}{\mathit{M}}_{\text{BC}}^0={\mathit{M}}_{\text{BC}}& {\mathit{C}}_{\text{BC}}^0={\mathit{C}}_{\text{BC}}& {\mathit{K}}_{\text{BC}}^0={\mathit{K}}_{\text{BC}}\\ {\mathit{M}}_{\text{CB}}^0={\mathit{M}}_{\text{CB}}& {\mathit{C}}_{\text{CB}}^0={\mathit{C}}_{\text{CB}}& {\mathit{K}}_{\text{CB}}^0={\mathit{K}}_{\text{CB}}\\ {\mathit{M}}_{\text{BB}}^0={\mathit{M}}_{\text{BB}}& {\mathit{C}}_{\text{BB}}^0={\mathit{C}}_{\text{BB}}& {\mathit{K}}_{\text{BB}}^0={\mathit{K}}_{\text{BB}}\end{array}\}.

When the free-field is known, the seismic equivalent loads can be obtained according to Eqs. (3) and (4), as shown in Eq. (5),

(5)\begin{array}{rl}& \left\{\begin{array}{c}0\\ 0\\ {\mathit{F}}_{\text{B}}\end{array}\right\}=\left[\begin{array}{ccc}{\mathit{M}}_{\text{II}}^0& {\mathit{M}}_{\text{IC}}^0& 0\\ {\mathit{M}}_{\text{CI}}^0& {\mathit{M}}_{\text{CC}}^0& {\mathit{M}}_{\text{CB}}\\ 0& {\mathit{M}}_{\text{BC}}& {\mathit{M}}_{\text{BB}}\end{array}\right]\left\{\begin{array}{l}{\ddot{\mathit{u}}}_{\text{I}}^0\\ {\ddot{\mathit{u}}}_{\text{C}}^0\\ {\ddot{\mathit{u}}}_{\text{B}}^0\end{array}\right\}\\ +& \left[\begin{array}{ccc}{\mathit{C}}_{\text{II}}^0& {\mathit{C}}_{\text{IC}}^0& 0\\ {\mathit{C}}_{\text{CI}}^0& {\mathit{C}}_{\text{CC}}^0& {\mathit{C}}_{\text{CB}}\\ 0& {\mathit{C}}_{\text{BC}}& {\mathit{C}}_{\text{BB}}\end{array}\right]\left\{\begin{array}{l}{\dot{\mathit{u}}}_{\text{I}}^0\\ {\dot{\mathit{u}}}_{\text{C}}^0\\ {\dot{\mathit{u}}}_{\text{B}}^0\end{array}\right\}+\left[\begin{array}{ccc}{\mathit{K}}_{\text{I}}^0& {\mathit{K}}_{\text{IC}}^0& 0\\ {\mathit{K}}_{\text{Cl}}^0& {\mathit{K}}_{\text{CC}}^0& {\mathit{K}}_{\text{CB}}\\ 0& {\mathit{K}}_{\text{BC}}& {\mathit{K}}_{\text{BB}}\end{array}\right]\left\{\begin{array}{l}{\mathit{u}}_{\text{I}}^0\\ {\mathit{u}}_{\text{C}}^0\\ {\mathit{u}}_{\text{B}}^0\end{array}\right\}.\end{array}

The equivalent input seismic load acting on the artificial boundary can be obtained from Eq. (5), as shown in Eq. (6),

(6){\mathit{F}}_{\text{B}}={\mathit{M}}_{\text{BC}}{\ddot{\mathit{u}}}_{\text{C}}^0+{\mathit{M}}_{\text{BB}}{\ddot{\mathit{u}}}_{\text{B}}^0+{\mathit{C}}_{\text{BC}}{\dot{u}}_{\text{C}}^0+{\mathit{C}}_{\text{BB}}{\dot{\mathit{u}}}_{\text{B}}^0+{\mathit{K}}_{\text{BC}}{\mathit{u}}_{\text{C}}^0+{\mathit{K}}_{\text{BB}}{\mathit{u}}_{\mathbf{B}}^0 .

The specific derivation process can be seen in Refs. [39,40]. It is evident that the equivalent seismic load, F_B, depends only on the motion of nodes on the substructure of artificial boundary and has no relationship with the motion of the remaining internal nodes. Therefore, after the motion of nodes on the substructure of artificial boundary is obtained by the FK method, the equivalent seismic load, F_B, is obtained by the artificial boundary substructure method (i.e., Eq. (6)). Then, FEM is used to establish the local domain model. After F_B is applied to the artificial boundary nodes of local FEM model, the simulation of local domain model could be performed to assess the seismic behavior of underground structure.

2.3 Step-by-step procedures of the hybrid approach

The FK can solve the free-field response under the action of dislocation point source. The response of corresponding nodes can be obtained with the source parameters, as well as velocity structure and coordinates of required nodes as the global model input. Since the finite element software ABAQUS has a friendly graphical user interface and can be processed in batch with Python, ABAQUS is selected to simulate the local domain model. The FK program, ABAQUS, and Python are implemented during the entire process of the FK-FEM hybrid approach. Note that any FEM software can be used as long as it provides the functionality to simultaneously add external loads to a number of nodes. The specific implementation process is as follows, and the step-by-step procedures is shown in Fig.2.

Step 1: Establish the global domain model and the local domain model. The calculation model is divided into the global domain and the local domain, and the global domain model and local domain model are built by FK program and ABAQUS, respectively.

Step 2: Establish the artificial boundary substructure model and export all nodes coordinates. Using the element birth and death of ABAQUS, internal elements of the local domain model are killed and the artificial boundary substructure model is obtained. The coordinates of all nodes on the substructure of artificial boundary are exported using Python script.

Step 3: Obtain all nodes displacements of the substructure of artificial boundary using FK program. The coordinates of all nodes obtained in Step 2 are fed into the FK program for calculation and all nodes’ displacements on the substructure of artificial boundary could be extracted.

Step 4: Obtain the equivalent seismic load F_B. The boundary of the substructure of artificial boundary is set as 3D viscous-spring artificial boundary [43] using Python script and the displacements obtained in Step 3 are applied to the corresponding nodes, respectively. ABAQUS is used to obtain the reaction force on artificial boundary nodes, that is, the equivalent seismic load, F_B, and export, F_B, using Python script.

Step 5: Simulate the dynamic response of the local region model for seismic evaluation of structures. The boundary of the local domain is set as viscous-spring artificial boundary and the equivalent seismic load, F_B, obtained in Step 4 is applied to the corresponding artificial boundary nodes. The process of FK-FEM hybrid approach is completed.

3 Verification of the FK-FEM hybrid approach

In this section, the calculation results by FK method and FK-FEM hybrid approach are compared to verify the correctness of the FK-FEM hybrid approach.

3.1 Verification of free-field

To verify the correctness of the FK-FEM hybrid approach in calculating free-field, the FK method and FK-FEM hybrid approach are employed to simulate free-field site response in the same homogeneous half-space. The model is shown in Fig.3, where the medium density ρ = 1800 kg/m³, P-wave velocity V_P = 1400 m/s, S-wave velocity V_S = 600 m/s, and damping ratio ν = 0.00625. The size of local FEM domain is 1000 m × 500 m × 500 m and its grid size is 20 m. The point source is located directly below the FEM model and its buried depth is 10 km. The seismic moment M₀ = 10¹³ N·m, i.e., moment magnitude M_W = 2.6. The relationship between M₀ and M_W is shown in Eq. (7). The strike, dip, and slip angle are 135°, 70°, and 30°, respectively. The source time function S(t) is a bell function, as shown in Eq. (8). Fig.4 shows the displacement results at points A (0.0, 0.0, 0.0), B (0.0, 0.0, 0.2), and C (0.0, 0.0, 0.4) from the FK method and FK-FEM hybrid approach, respectively. It is evident that the simulation results from the two methods are in good agreement, which indicate correctness of FK-FEM hybrid approach.

(7)M_{\text{W}} =\frac{2}{3}\lg (M_{\text{0}} -9.1) ,

(8)

3.2 Verification of scattering-field

To verify the correctness of FK-FEM hybrid approach in calculating scattering-field, a block embedded in a homogeneous half-space is set with the same parameters as the half-space. The displacements of 3 points in the block are calculated by FK-FEM hybrid approach and compared with the displacements by FK method. The model is illustrated in Fig.5. The block size is 20 m × 12 m × 12 m and the buried depth is 10 m. Tie constraint is applied to the contact between block and half-space. The homogeneous half-space parameters are: density ρ = 3000 kg/m³, P-wave velocity V_P = 2400 m/s, elastic modulus E = 11.52 × 10³ MPa, and Poisson’s ratio v = 0.3333. The incident wave is a P wave with incident angle 30° and the waveform is Ricker wave. Fig.6 shows the displacements at points A, B, and C in this model. The results by the two methods are in good agreement, which verifies the correctness of FK-FEM hybrid approach.

4 Simulation examples and results

In this section, the FK method is employed to evaluate the displacements response of 8 adjacent points in three directions, which further shows that this method can effectively capture the spatial variability characteristics of seismic motion input. Then, taking Daikai subway station as an example, the seismic response of an underground structure from a point dislocation source is simulated using FK-FEM hybrid approach, after which the effects of crustal velocity structure and source parameters on the seismic response of Daikai station are discussed.

4.1 Spatial variability of seismic motion

To illustrate that the spatial variability characteristics of seismic motion can be captured by this approach, a layered half-space is built, as shown in Fig.7. The crustal velocity structure is shown in Tab.1. A point source is located 20 km deep from the surface. The seismic moment M₀ = 4.6 × 10¹⁹ N·m, i.e. M_W = 7.0, and the strike (ϕ), dip (δ), and slip (λ) angle are 0°, 45°, and 45°, respectively. Then, the moment tensor can be obtained by referring to Eq. (9). The slip rate function is a bell-shaped function, as shown in Eq. (8). Fig.7 shows the displacement time history at 8 vertices of the cuboid. The size of the site in this example is 200 m × 200 m × 100 m, which generally satisfies the seismic analysis of common underground structures.

(9)\{\begin{array}{rl}& M_{xx}=-M_0(\sin \delta \cos \lambda \sin 2\varphi +\sin 2\delta \sin \lambda {\sin }^2\varphi ),\\ & M_{yy}=M_0(\sin \delta \cos \lambda \sin 2\varphi -\sin 2\delta \sin \lambda {\cos }^2\varphi ),\\ & M_{\mathit{\text{zz}}}=M_0(\sin 2\delta \sin \lambda ),\\ & M_{xy}=M_0\left(\sin \delta \cos \lambda \cos 2\varphi +\frac{1}{2}\sin 2\delta \sin \lambda \sin 2\varphi \right),\\ & M_{x\mathit{\text{z}}}=-M_0(\cos \delta \cos \lambda \cos \varphi +\cos 2\delta \sin \lambda \sin \varphi ),\\ & M_{y\mathit{\text{z}}}=-M_0(\cos \delta \cos \lambda \sin \varphi -\cos 2\delta \sin \lambda \sin \varphi ).\end{array}

By observing Fig.7, pronounced differences can be seen between the displacement amplitude and the direction of the maximum displacement at different calculation points. For example, in the y-direction, the maximum displacements at vertices A, D, E, and F are on the positive side, but the maximum displacements at other vertices are negative. In the x-direction, the maximum displacement of points B and H are 1.106 and 1.645 m, respectively, and the amplitude difference is 48.73%. The results further verify that the input seismic wave field cannot be simply assumed as plane wave in the seismic response analysis of underground structures. In addition, it shows that the FK-FEM hybrid approach can capture the spatial variability of seismic motion caused by source and propagation path.

4.2 Influence of crustal velocity structure on underground structure

Reflection and refraction occur in the process of seismic wave propagation in crustal layers, resulting in changes in the amplitude and direction of seismic wave.

To explore the influence of crustal velocity structure on the seismic response of underground structures, three crustal velocity structure cases with the same surface soil parameters are established, as shown in Tab.2. Compared to Case 2, the crust of Case 1 is softer, while the crust of Case 3 is harder. The seismic moment M₀ = 4.6 × 10¹⁹ N·m, i.e., M_W = 7.0. The point source to the center of the local domain model is 40 km and the buried depth of the point source is 20 km, as shown in Fig.8. The strike (ϕ), dip (δ), and slip (λ) angle are 52°, 88°, and 176°, respectively. The slip rate function is Eq. (8).

Taking Daikai subway station as the prototype, ABAQUS is used to build a 3D local domain model. The buried depth of the subway station is 5 m, and the finite element model size and Subsection size of station are shown in Fig.8. 3D viscous-spring artificial boundary are employed as the boundary conditions. Due to the complexity of the 3D subway station, elastic material model is employed for concrete to improve the computational efficiency. The elastic modulus of concrete is 3 × 10⁴ MPa, density is 2600 kg/m³, and Poisson’s ratio is 0.2. Solid hexahedron element (C3D8) is used to model the soil and subway station. The grid sizes of soil and station are 5 m × 5 m × 5 m and 1 m × 1 m × 1 m, respectively, and the grid size of columns inside the station is 0.2 m × 0.2 m × 0.2 m.

Tie constraint is adopted for the contact between soil and structure. As the station has many cross sections, this paper only takes Cross Sections 1, 2, and 3 as examples for analysis.

Fig.9 shows the maximum relative displacements of the side-wall at the three sections (i.e., the maximum displacement of all nodes on the side wall relative to the bottom of the side wall). Fig.10 shows the time histories of relative interlayer displacement angle, as well as the deformation diagram when the relative interlayer displacement angle is the largest. It should be noted that only results in the y-direction are presented in Fig.9 and Fig.10 since the deformation of the structure in the y-direction is the largest.

According to Fig.9, the maximum relative displacement of the side wall at Cross Section 1 is different from those of Cross Sections 2 and 3. However, the maximum deformation for the three crustal velocity structures at each cross section follows the same pattern. This may be due to the fact that the deformation characteristics are determined by the stiffness characteristics of the structure. In addition, it is evident that the softer the crustal velocity structure, the greater the maximum deformation of the structure. This is especially evident at Cross Section 3, where the maximum displacement of the side wall of Case 1 is 117% larger than that of Case 3.

Fig.10 illustrates that the softer the crustal velocity structure, the greater the maximum interlayer displacement angle, which follows the same pattern as results in Fig.9. This was also evident at Cross Section 3, where the maximum displacement angle of the side wall of Case 1 is 96.5% larger than that of Case 3. In addition, it is observed that with the softening of the crustal velocity structure, the time at which the maximum interlayer displacement angle occurs is obviously delayed. For example, at Cross Section 1, the time at which the maximum interlayer displacement angle occurs under the crustal velocity structure Case 1 is delayed 3.25 s compared to that of Case 3.

Therefore, crustal velocity structure will significantly change the deformation level and time at which the maximum deformation of the structure occurs. The softer the crustal velocity structure, the more intense the underground structure deformation and more delayed the maximum deformation. The influence of the crustal velocity structure should also be considered in the seismic design of underground structures.

4.3 Influence of source parameters on underground structure

The source parameters directly determine the characteristics of seismic wave field, resulting in differences in the seismic response of underground structures.

To further illustrate the impact of source parameters on underground structures, three groups of strike, dip, and slip angles are set based on the model in Subsection 4.2 and with reference to Ref. [44], as shown in Tab.3. The other source parameters are the same as those in Subsection 4.2. Fig.11 and Fig.12 show the time histories of relative interlayer displacement angle in x- and y-directions at Cross Sections 1, 2, and 3, respectively.

Fig.11 and Fig.12 show that the maximum interlayer displacement angle of the structure changes noticeably under different source parameters. In general, the maximum interlayer displacement angles in the y-direction are greater than those in the x-direction. However, the maximum interlayer displacement angles in the x-direction at the top layer of Cross Section 3 of Source 2 and Source 3 are greater than that in the y-direction.

In addition, at the top layer of Cross Section 3, the maximum interlayer displacement angles in the x-direction of Source 1 and Source 3 are negative (–0.2220‰ and –0.6178‰, respectively), but the maximum interlayer displacement angle of Source 2, is positive (0.3272‰). The time at which the maximum interlayer displacement angle occurs under the Source 2 is evidently later than that of Source 1 and Source 3, which indicates that the size and direction of the maximum deformation have changed.

Therefore, source parameters significantly affect the maximum relative interlayer displacement angle, direction of structural damage, and time at which the maximum deformation of the structure occurs.

5 Conclusions

In this paper, the FK-FEM hybrid approach is established based on the domain reduction idea and artificial boundary substructure method. FK is employed to simulate the seismic wave propagation from source, along wave propagation path to ground surface. FEM is used to simulate the seismic responses of local sites and underground structures. Therefore, the physics-based simulations from point dislocation source to underground structure could be performed with the FK-FEM hybrid approach. By comparing with the results of FK, the correctness of the FK-FEM approach is verified. The main conclusions are as follows.

1) The FK-FEM hybrid approach is different from the plane wave input method and its seismic wave field input contains information, such as source and propagation path, as well as spatial variability characteristics caused by propagation path.

2) The velocity structure of crust has significant influence on the deformation of the underground structure and time at which the maximum deformation occurs. The softer the velocity structure of crust, the more unfavorable to underground structure deformation and more delayed the maximum deformation. In the example of this paper, the influence of velocity structure on the maximum interlayer displacement angle of the underground structure is 96.4%. The influence of crustal velocity structure on underground structure should also be considered in seismic design.

3) Source parameters not only significantly affect the maximum deformation of the underground structure and time at which the maximum deformation occurs, but also change the direction of structural damage. Further research is required to clarify the specific influence mechanism.

4) The simulation in this paper only takes into consideration an example of an underground structure. It is worth noting that FK-FEM hybrid approach can also be used for aboveground structures and an appropriate constitutive model can be used in the local area model to consider the nonlinearity of surface soil and structure.

References

[1] Dong Y H, Peng F L, Guo T F. Quantitative assessment method on urban vitality of metro-led underground space based on multi-source data: A case study of Shanghai Inner Ring area. Tunnelling and Underground Space Technology, 2021, 116: 104108

[2] Xie H, Zhang Y, Chen Y, Peng Q, Liao Z, Zhu J. A case study of development and utilization of urban underground space in Shenzhen and the Guangdong-Hong Kong-Macao Greater Bay Area. Tunnelling and Underground Space Technology, 2021, 107: 103651

[3] Guo D, Zhu X, Xie J, Zhang C, Zhao Z, Zheng L. Undergraduate program for urban underground space engineering in China: Exploration and practice. Tunnelling and Underground Space Technology, 2021, 116: 104084

[4] Zaini F, Hussin K, Raid M M. Legal considerations for urban underground space development in Malaysia. Underground Space, 2017, 2(4): 234–245

[5] Yang J, Zhuang H, Wang W, Zhou Z, Chen G. Seismic performance and effective isolation of a large multilayered underground subway station. Soil Dynamics and Earthquake Engineering, 2021, 142: 106560

[6] Iida H, Hiroto T, Yoshida N, Iwafuji M. Damage to Daikai subway station. Soils and Foundations, 1996, 36(Special): 283–300

[7] Samata S, Ohuchi H, Matsuda T. A study of the damage of subway structures during the 1995 Hanshin-Awaji earthquake. Cement and Concrete Composites, 1997, 19(3): 223–239

[8] Yamaguchi A, Mori T, Kazama M, Yoshida N. Liquefaction in Tohoku district during the 2011 off the Pacific Coast of Tohoku Earthquake. Soil and Foundation, 2012, 52(5): 811–829

[9] Xu Z, Du X, Xu C, Hao H, Bi K, Jiang J. Numerical research on seismic response characteristics of shallow buried rectangular underground structure. Soil Dynamics and Earthquake Engineering, 2019, 116: 242–252

[10] Li W, Chen Q. Effect of vertical ground motions and overburden depth on the seismic responses of large underground structures. Engineering Structures, 2020, 205: 110073

[11] Hashash Y M A, Hook J J, Schmidt B, I-Chiang Yao J. Seismic design and analysis of underground structures. Tunnelling and Underground Space Technology, 2001, 16(4): 247–293

[12] Kuesel T R. Earthquake design criteria for subways. Journal of the Structural Division, 1969, 95(6): 1213–1231

[13] Liu J, Wang W, Dasgupta G. Pushover analysis of underground structures: method and application. Science China. Technological Sciences, 2014, 57(2): 423–437

[14] Penzien J. Seismically induced racking of tunnel linings. Earthquake Engineering & Structural Dynamics, 2000, 29(5): 683–691

[15] TateishiA. A study on loading method of seismic deformation method. Proceedings of the Civil Society, 1992, 1992(441): 157–166 (in Japanese)

[16] Tateishi A. A study on seismic analysis methods in the cross section of underground structures using static finite element method. Structural Engineering/Earthquake Engineering, 2005, 22(1): 41s–54s

[17] WangJ N. Seismic Design of Tunnels: A Simple State-of-the-Art Design Approach. New York: Brinckerhoff Parsons Inc., 1993

[18] Zhao M, Gao Z, Du X, Wang J, Zhong Z. Response spectrum method for seismic soil-structure interaction analysis of underground structure. Bulletin of Earthquake Engineering, 2019, 17(9): 5339–5363

[19] Xu C, Zhang Z, Li Y, Du X. Validation of a numerical model based on dynamic centrifuge tests and studies on the earthquake damage mechanism of underground frame structures. Tunnelling and Underground Space Technology, 2020, 104: 103538

[20] Zhu T, Wang R, Zhang J M. Evaluation of various seismic response analysis methods for underground structures in saturated sand. Tunnelling and Underground Space Technology, 2021, 110: 103803

[21] IwatateTKobayashi YKusuHRinK. Investigation and shaking table tests of subway structures of the Hyogoken-Nanbu earthquake. In: 12th World Conference on Earthquake Engineering. Auckland: New Zealand Society for Earthquake Engineering Upper Hutt, 2020

[22] Zhang Z, Li Y, Xu C, Du X, Dou P, Yan G Y. Study on seismic failure mechanism of shallow buried underground frame structures based on dynamic centrifuge tests. Soil Dynamics and Earthquake Engineering, 2021, 150: 106938

[23] Trifonov O V, Cherniy V P. A semi-analytical approach to a nonlinear stress–strain analysis of buried steel pipelines crossing active faults. Soil Dynamics and Earthquake Engineering, 2010, 30(11): 1298–1308

[24] Wang X, Chen J, Xiao M. Seismic responses of an underground powerhouse structure subjected to oblique incidence SV and P waves. Soil Dynamics and Earthquake Engineering, 2019, 119: 130–143

[25] Wang X, Xiong Q, Zhou H, Chen J, Xiao M. Three-dimensional (3D) dynamic finite element modeling of the effects of a geological fault on the seismic response of underground caverns. Tunnelling and Underground Space Technology, 2020, 96: 103210

[26] Jiang J, Nggar H M, Xu C, Zhong Z, Du X. Effect of ground motion characteristics on seismic fragility of subway station. Soil Dynamics and Earthquake Engineering, 2021, 143: 106618

[27] Zhu Z, Tang Y, Ba Z, Wang K, Gong W. Seismic analysis of high-speed railway irregular bridge–track system considering V-shaped canyon effect. Railway Engineering Science, 2022, 30(1): 57–70

[28] Zhang L, Wang J, Xu Y, He C, Zhang C. A procedure for 3D seismic simulation from rupture to structures by coupling SEM and FEM. Bulletin of the Seismological Society of America, 2020, 110(3): 1134–1148

[29] He C H, Wang J T, Zhang C H, Jin F. Simulation of broadband seismic ground motions at dam canyons by using a deterministic numerical approach. Soil Dynamics and Earthquake Engineering, 2015, 76: 136–144

[30] Mccallen D, Petersson N, 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

[31] Yu H T, Yang Y S, Yuan Y, Duan K P, Gu Q. A comparison between vibration and wave methods in seismic analysis of underground structures. China Earthquake Engineering Journal, 2019, 41(04): 845–852

[32] 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

[33] BielakJ. Reply to “Comment on ‘Domain reduction method for three-dimensional earthquake modeling in localized regions, Part I: Theory,’ by J. Bielak, K. Loukakis, Y. Hisada, and C. Yoshimura, and ‘Part II: Verification and applications,’ by C. Yoshimura, J. Bielak, Y. Hisada, and A. Fernandez,” by E. Faccioli, M. Vanini, R. Paolucci, and M. Stupazzini. Bulletin of the Seismological Society of America, 2005, 95(2): 770−773

[34] Yoshimura C. Domain reduction method for three-dimensional earthquake modeling in localized regions, Part II: Verification and applications. Bulletin of the Seismological Society of America, 2003, 93(2): 825–841

[35] Li W, Assimaki D. Simulating soil stiffness degradation in transient site response predictions. Soil Dynamics and Earthquake Engineering, 2010, 30(5): 299–309

[36] Crempien J G F, Archuleta R J. UCSB method for simulation of broadband ground motion from kinematic earthquake sources. Seismological Research Letters, 2015, 86(1): 61–67

[37] Brun M, De Martin F, Richart N. Hybrid asynchronous SEM/FEM co-simulation for seismic nonlinear analysis of concrete gravity dams. Computers & Structures, 2021, 245: 106459

[38] Casadei F, Gabellini E, Fotia G, Maggio F, Quarteroni A. A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems. Computer Methods in Applied Mechanics and Engineering, 2002, 191(45): 5119–5148

[39] Liu J, Tan H, Bao X, Wang D, Li S. Seismic wave input method for three-dimensional soil-structure dynamic interaction analysis based on the substructure of artificial boundaries. Earthquake Engineering and Engineering Vibration, 2019, 18(4): 747–758

[40] Liu J, Bao X, Wang D, Wang P. Seismic response analysis of the reef-seawater system under incident SV wave. Ocean Engineering, 2019, 180: 199–210

[41] AkiKRichardsP G. Quantitative Seismology. 2nd ed. Sausalito: University Science Books, 2002

[42] Ba Z, Sang Q, Wu M, Liang J. The revised direct stiffness matrix method for seismogram synthesis due to dislocations: From crustal to geotechnical scale. Geophysical Journal International, 2021, 227(1): 717–734

[43] 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

[44] Fukuyama E, Ishida M, Dreger D, Kawai H. Automated seismic moment tensor determination by using on-line broadband seismic waveforms. Journal of the Seismological Society of Japan, 1998, 51(1): 149–156

RIGHTS & PERMISSIONS

Higher Education Press