#### Abstract

A new internal substructure method for seismic wave input in soil-structure systems was recently proposed. This method simplifies the calculation of equivalent input seismic loads and avoids the participation of artificial boundaries in the process of seismic wave input. However, in previous research and applications, the internal substructures are usually intercepted down from the free surface, which forms large substructures and increases the computational effort for data management on the substructure nodes, especially for deep underground structures. In this study, the internal substructure method is modified by intercepting the internal substructures entirely beneath the free surface and adjacently around the underground structures. Then, the equivalent input seismic loads are obtained through the dynamic analysis of the internal substructures and applied to the corresponding positions of the total soil-structure models. Thus, the earthquake energy can be more efficiently input into the region near the underground structures without losing computational accuracy. We provide the detailed implementation procedures of this modified method and validate its applicability and accuracy through the scattered problems of underground cavities in homogeneous and layered half-space sites.

#### 1. Introduction

With the rapid development of underground engineering in recent years, the seismic safety of underground structures has received increasing attention. The current seismic analysis methods for underground structures can be generally classified into experimental research [1, 2], theoretical study [3–5], and computational simulation [6–8]. Among them, computational analysis has the advantage of handling complex geometric shapes and material conditions; thus, it becomes the most practical approach for the seismic analysis of underground structures.

To simulate the seismic response of underground structures, which constitutes a semi-infinite system with the unbounded foundation, we commonly establish a finite near-field model including the underground structures and surrounding soil and set the artificial boundaries on the cutoff boundaries to simulate the wave radiation into infinity. The representative achievements in the artificial boundaries include the boundary element method [9, 10], the infinite element method [11], perfectly matched layers [12, 13], transmitting boundaries [14, 15], viscous boundaries [16], viscoelastic boundaries [17, 18], and exact absorbing boundaries [19, 20]. However, due to the application of the artificial boundaries, a new problem appears in inputting the seismic waves into the finite model without affecting the transmission of outgoing waves. To settle this issue, the commonly used schemes, represented by the domain reduction method (DRM) [21–24] and the boundary reaction method (BRM) [25], generally convert the incident waves into the equivalent input seismic loads. DRM is a time domain method which can simulate the earthquake source and propagation path effect. And, BRM is a hybrid frequency-time domain two-step method which separates the simulated system into the far-field linear domain and the near-field domain that can consider the nonlinear mechanical behavior. In the 1990s, our research team proposed the wave method [26] based on the viscoelastic artificial boundaries. This method has been demonstrated to be precise in a multitude of applications [6, 27–29], but it is complicated in calculation. Recently, based on the wave method concept, we proposed the boundary substructure method (BSM) [30] for seismic wave input. This method only needs to calculate the equivalent input seismic loads on one layer of nodes through a dynamic analysis of the artificial boundary substructures, so it is more simplified in implementation. However, both methods have a weakness: the participation of artificial boundaries in the process of seismic wave input. On account of this feature, only the stress-type artificial boundaries can be applied. In addition, since the cutoff boundaries are usually set far from the underground structures to reduce the approximation error for wave absorption of artificial boundaries, the numerical models of soil-structure systems are generally exceedingly large. The equivalent input seismic loads should also be separately calculated on each boundary node, which leads to a heavy computational burden. To overcome these defects, we recently extended BSM to more general situations and proposed a more universally applicable method, namely, the internal substructure method (ISM) [31]. In this method, the substructure model is directly intercepted from the original soil-structure model, which is of smaller size and higher computational efficiency. In addition, this method avoids the participation of the artificial boundaries in the process of inputting seismic waves. Therefore, the restrictions of the traditional methods on the type of artificial boundaries are completely removed.

However, in the previous research and applications of ISM, the internal substructures are usually intercepted down from the free surface [31]. For shallow underground structures, this scheme is acceptable, since the structures are near the surface. For deep underground structures, unsurprisingly, this method leads to a large internal substructure, which in turn increases the effort for data management on the substructure nodes. In this study, we modify the ISM by defining the region of interest completely below the free surface and intercepting the enclosed substructures around the underground structures. Thus, the seismic waves can be effectively input near the deep underground structures, and the computational workloads can be significantly reduced. The implementation procedures of the modified method are introduced, and the scattered problems of underground cavities in the homogeneous and layered half-space sites are used as numerical cases to verify the effectiveness of this modified internal substructure method.

#### 2. Modified Internal Substructure Method

##### 2.1. Theory

Similar to the BSM and original ISM, to input the seismic waves into the numerical model, we should convert the wave motions into the equivalent input seismic loads. Here, in the modified ISM, we define the region of interest, where the soil and underground structures interact with each other, completely beneath the free surface (see the region in the blue node box in Figure 1(a)). The seismic waves should be input around the region of interest. According to the concept of the ISM, the correct equivalent input seismic loads should ensure that the motions at the input locations are identical with those of the free wave field. Therefore, the equivalent input seismic loads can be obtained through the dynamic analysis of the corresponding free-field model, as shown in Figure 1(b).

**(a)**

**(b)**

**(c)**

We classify the nodes in the finite element model into five categories according to their positions: interior nodes (marked by I), exterior nodes (marked by E), and three layers of internal substructure nodes (marked by A, B, and C from the outside to the inside). Notably, to realize the consistency of the submatrixes and the correspondence of the seismic wave input locations, the mesh generation connected to nodes E, A, B, and C should be identical in the soil-structure model and free-field model.

According to the above-mentioned nodal classification, the motion equation of the free-field finite element model can be obtained as follows:where **K**, **C**, and **M** are the stiffness, damping, and mass matrixes, respectively; **F** is the vector of nodal force; **u**, , and are the vectors of the displacement, velocity, and acceleration, respectively; the subscripts E, A, B, C, and I denote the aforementioned nodal classification.

We divide the free-field model into the interior and exterior regions by nodes A. According to the theory of the isolator, the wave motion inside the interior region can be determined uniquely once the nodal forces of exterior region are identified. Here, the wave field inside the interior region is supposed to be the free wave field, namely,where superscript 0 denotes the free-field motion.

For the exterior region, since the input seismic waves are neglected in this domain, nodes A and E remain stationary during the simulation, i.e.,

By substituting equation (3) into equation (1), we obtain the equivalent input loads as follows:

The values of **F**_{A} and **F**_{B} only involve the free-field motions of nodes B and C and the submatrixes with subscripts AB, BB, and BC. In order to reduce the calculation region, we establish a finite element model of internal substructure, as shown in Figure 1(c), which is composed of nodes A, B, and C, as well as the two layers of elements between them. The motion equation of the internal substructure model can be obtained:where the superscript denotes the internal substructure model.

According to the theory of finite element discretization, if the grid partition and the material properties of the internal substructure model are identical with those of the free-field model, we can obtain the following equations:

After substituting equation (6) into equation (5), we fix the A nodes (), apply the free-field motions on the B and C nodes ( and ), and then obtain the following equations:

By comparing equation (4) and equation (7), it can be found that the nodal forces on the A and B nodes calculated by the free-field model and the internal substructure model are the same. Therefore, the equivalent input seismic loads can be gained through the dynamic analysis of the internal substructure model under the imposed free-wave-field motions.

It is important to emphasize that the derivation the proposed method is based on the finite element method (FEM). Therefore, the accuracy of this method is identical with that of FEM. According to the previous research [32], for the wave propagation simulation problems, the mesh size should meet the following condition:where *c*_{s} is the small-strain shear-wave velocity in the media and *f*_{cut} is the cutoff frequency.

Normally, in the practice of wave propagation simulation, equation (8) should be checked before calculation. However, if the seismic wave leads to high shear strains and significant softening of the soil, meshes may be required to be nicer [33].

##### 2.2. Implementation Procedures

Based on the aforementioned theoretical derivations, the original ISM can be modified to optimize the seismic wave input of the underground structure-soil systems. Then, we introduce the implementation procedures of this modified ISM as follows:(1)Establish the soil-structure finite element model, as shown in Figure 1(a), and determine the spatial range of seismic wave inputs according to the position of the underground structures. To reduce the computational effort, this spatial range should be adjacently around the underground structures and entirely beneath the free surface. Then, we can establish the internal substructure model, as shown in Figure 2(a), by deleting nodes E and I of the original soil-structure model and the finite elements connected to these nodes.(2)For any selected incident waves, the corresponding free wave motions can be obtained through the free-field analysis algorithms [34, 35]. Then, by applying the free-wave motions **u**_{B} and **u**_{C} on the corresponding nodes B and C of the internal substructure model and fixing nodes A, as shown in Figure 2(a), we can conduct a dynamic analysis on the internal substructure. The reaction forces on nodes A and B can be gained, which are the equivalent input seismic loads **F**_{A} and **F**_{B}.(3)Apply the equivalent input seismic loads **F**_{A} and **F**_{B} on the corresponding nodes A and B of the original soil-structure model, as shown in Figure 2(b), and conduct a dynamic analysis on this soil-structure model. Hence, the seismic waves are effectively input into the soil-structure model.

**(a)**

**(b)**

##### 2.3. Characteristics

In this section, we will compare this modified method with original ISM, BSM, and the famous DRM. In this way, the characteristics of the modified ISM can be concluded and its scope of application can be defined.

In this modified ISM, the enclosed internal substructures, which are entirely beneath the free surface, are intercepted adjacently around the underground structures. Thereby, the size of the substructure model becomes significantly smaller than that of the original ISM, where the substructure model is intercepted down from the free surface. This size reduction feature effectively reduces the computational effort, especially the tremendous data management on the input locations.

In BSM, the artificial boundaries are included in the substructure. Therefore they participate in the calculation of equivalent seismic loads, and thus the restrictions on the applied artificial boundaries are more severe in BSM. Namely, they need to be the stress-type boundary conditions (such as the viscous boundaries [16] and the viscoelastic boundaries [17, 18]), while the displacement-type artificial boundaries such as transmitting boundaries [14, 15] or the perfectly matched layers [12, 13] are not applicable in BSM. The original and modified ISMs avoid the participation of the artificial boundaries in the process of equivalent loads calculation by spatially separating the input of seismic waves and the absorption of scattered waves. In this way, the limitations on the type of artificial boundaries are completely released. Therefore, ISM is more universal in application.

However, it is necessary to note that although the sizes of substructures get much smaller in the modified and original ISM, we should calculate the equivalent seismic loads on two layers of nodes. When the calculation domain shows strong nonlinearity or we are interested in the overall results of the calculation domain, the entire wave field needs to be obtained, rather than only the seismic response of the structures and the scattered wave field in the external region. In this situation, the seismic waves should be input through the artificial boundaries by BSM, on account of its advantage in reducing the data management from two-layer nodes to only one-layer nodes.

Now, we should discuss the relationship between the famous DRM and ISM. The differences between these two methods can be found in the theoretical principles, the mechanical models, and the calculation times, as well as the locations of seismic wave input:(1)DRM and ISM involve the conversion of the stiffness matrixes related to the external region and the internal region, respectively. In DRM, the submatrixes related to the external region are converted to submatrixes related to the interface between the internal region and the external region, according to the force equilibrium in the interface. While in the original and modified ISM, the nature of FEM is applied, namely, the nodal forces are only transmitted through the connected elements. If there is no element between two nodes, there is no interaction force. By introducing nodes C, the submatrixes related to the internal region are converted into the submatrixes between nodes B and C. In this way, the internal submatrixes are eliminated in the formulas of the equivalent seismic loads.(2)In the practical applications of DRM, in order to obtain the equivalent seismic loads, it is necessary to conduct two dynamic calculations under different boundary conditions: for , we should apply the free-wave-field motions on the nodes “e” and fix the nodes “b”. As for , the free-wave motions should be applied on the nodes “b” and the nodes “e” should be fixed. While in the original and modified ISM, only an internal substructure model, as shown in Figure 1c, is established. And, the equivalent seismic loads can be obtained through only one-time dynamic calculation of the internal substructure, under the free-wave-field motions applied on nodes B and C.(3)Since there are no limitations on the properties of the elements between nodes A and B in ISM, when the substructure is extended to the artificial boundaries, the properties of the elements between nodes A and B can be replaced by the properties of the artificial boundaries, and ISM degenerates into BSM. In other words, ISM can accomplish the seismic wave input on any location of the calculation system, including the artificial boundaries. However, DRM cannot deal with the seismic wave input on the artificial boundaries.

#### 3. Verification and Application of the Modified ISM

##### 3.1. Underground Cavity in the Homogeneous Half-Space Site

In this section, we validate the applicability of the modified ISM in the seismic analysis of underground structures through several numerical cases. First, we consider the wave scattering problem of the underground cavity in a homogeneous half-space site. The computational model is shown in Figure 3. The model dimensions are 200 m × 100 m. The circular cavity locates in the center of the model, and its diameter is 6 m. *θ* is the incident angle of the seismic wave. The material parameters of half space are as follows: mass density *ρ* = 2000 kg/m^{3}, shear-wave velocity *c*_{S} = 200 m/s, and Poisson’s ratio . The four-node plane strain elements are adopted in this model, and the consistent viscoelastic artificial boundary elements [36] are set on the cutoff boundaries. The red and blue dashed boxes in Figure 3 denote the two internal substructures intercepted according to the modified and original ISMs, which are denoted as models 1 and 2, respectively. The dimensions of these two internal substructures are 50 m × 50 m and 50 m × 75 m. Meanwhile, the accuracy of the modified method is verified by comparing the results with those of BSM, which is denoted as model 3. Observation point A is the midpoint on the free surface, and points B and C are the top and bottom points of the cavity, respectively. A plane SV impulse wave with the pulse duration of 0.2 s is normally input into the model, whose displacement time history is plotted in Figure 4.

Figure 5(a) shows the displacement on point A calculated by the modified ISM through model 1. The summation of this displacement time history and free-field motion on the same point, as shown in Figure 5(b), is consistent with the results calculated by the original ISM and BSM. Because the total wave field can be divided into the scattered wave field and free wave field, the above comparison can demonstrate that the wave field outside the internal substructures calculated by the modified ISM is exactly the scattered wave field. Figures 5(c) and 5(d) show the displacements on nodes B and C. In the internal substructures, these three numerical methods have consistent results. Therefore, the total wave field in the wave input region can be precisely calculated by the modified ISM.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 6 shows the snapshots of the displacement wave field. The dashed boxes denote the internal substructures; the red dashed box is adopted to input the seismic waves. It is important to emphasize that the wave field distributions in the red dashed box calculated by ISM are consistent with those of BSM, which indicates that the present method can efficiently input the seismic waves into the internal region. In addition, there are only the scattered waves in the region outside the internal substructures, and these scattered wave fields are consistent with those of the BSM (model 3).

**(a)**

**(b)**

**(c)**

##### 3.2. Underground Cavity in the Layered Half-Space Site

The layered half-space sites and oblique incidence of the seismic wave are more practical in engineering cases. Here, we consider a two-layer site with an underground cavity, as shown in Figure 7. The dimension of the model is 200 m × 200 m, and the cladding is 100 m thick. The cavity with a diameter of 6 m locates in the middle of the material interface. The material parameters of the elastic half space and cladding are mass density *ρ*_{1} = *ρ*_{2} = 2000 kg/m^{3}, velocity of shear wave *c*_{S1} = 200 m/s, *c*_{S2} = 100 m/s, and Poisson’s ratio . The dimensions of the internal substructures of model 1 and model 2 are 80 m × 80 m and 80 m × 140 m, respectively. Similar to the previous case, the boundary substructure method is adopted and denoted as model 3. Observation points A and B are the top and bottom point of the cavity, respectively. The pulse plane SV wave in Figure 4 is obliquely input into the model, and the incident angle is *θ* = 20°.

The displacements on nodes A and B are plotted in Figure 8. The results calculated by different methods are consistent with each other. Figure 9 shows the snapshots of the displacement wave fields. As expected, the wave field distributions in the internal regions (red dashed boxes in Figure 9) calculated by the ISM are consistent with those of the BSM.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 3.3. Deep Underground Cavity in the Layered Half-Space Site

To further demonstrate the superiority of the modified ISM for the seismic wave input in deep underground structure-soil systems, in this third numerical example, we consider a deep buried circular cavity in a two-layer site, as shown in Figure 10. The dimension of the model is 200 m × 600 m, and the cladding is 100 m thick. The diameter of the cavity is 30 m, and its depth is 500 m. The material parameters of the model are mass density *ρ*_{1} = 2500 kg/m^{3}, *ρ*_{2} = 2200 kg/m^{3}, velocity of shear wave *c*_{S1} = 800 m/s, *c*_{S2} = 500 m/s, and Poisson’s ratio *ν*_{1} = *ν*_{2} = 0.25, where subscripts 1 and 2 denote the elastic half space and the cladding. The red dashed box, which is marked as model 1, denotes the internal substructure adopted in the modified ISM, which is entirely beneath the free surface, and its dimension is 40 m × 40 m. Model 2 is the substructure in the original ISM, which is intercepted down from the free surface, and its dimension is 40 m × 520 m. Similarly, the boundary substructure method is adopted and denoted as model 3. The depths of observation points A and B are 480 m and 515 m, respectively. The pulse plane SV wave in Figure 4 is vertically input into the model.

The displacements in the *x*-direction at nodes A and B are plotted in Figure 11. Obviously, the results calculated by different methods are consistent with one another. Furthermore, we provide the peak displacement and peak von Mises stress at points A and B calculated by different methods in Table 1. The relative errors between the results of the ISMs and BSM are calculated and listed in the round brackets in Table 1. Compared to the BSM, the relative error of the modified ISM is less than 1%, which is even slightly more accurate than the original ISM.

**(a)**

**(b)**

Table 2 provides the quantities of elements and nodes in different substructure models. The model size significantly decreases through the modified ISM. Specifically, the numbers of elements and nodes decrease by 85.2% and 86.0% compared to the original ISM and by 88.6% and 89.2% compared to the BSM. This size reduction effectively decreases the computational workloads and further enhances the analysis efficiency.

Figure 12 shows the snapshots of the displacement wave fields. All wave field distributions in the internal regions (red dashed boxes in Figure 12) calculated by the internal substructure method are consistent with those of the BSM. Therefore, the accuracy and effectiveness of the modified method for the seismic wave input in deep underground structure-soil systems are demonstrated.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Conclusions

The internal substructure method has the advantages of reducing the model size and avoiding the participation of artificial boundaries in the seismic wave input. In the previous research and applications, however, the internal substructures are usually intercepted down from the free surface, which leads to relatively large substructures. The workloads for data management on the substructure nodes, especially for deep underground structures, also remain heavy. In this study, we modify the internal substructure method by defining the internal substructures completely beneath the free surface and intercepting these enclosed substructures adjacently around the underground structures. The equivalent input seismic loads can be obtained through the dynamic analysis of the internal substructures and applied to the total soil-structure models. Thus, the earthquake energy can be input into the region near the underground structures with less computational effort. The application of this modified method is studied through the wave scattering problems of underground cavities in homogeneous and layered half-space sites. The results demonstrate the high calculation precision of the modified ISM. Moreover, the superiority of the modified ISM in computational efficiency is prominent when the structures are deeply buried.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Acknowledgments

This study was supported by the National Natural Science Foundation of China (grant no. 51878384) and National Key Research and Development Program of China (grant no. 2018YFC1504305). Financial support from these organizations is gratefully acknowledged.