Abstract
Transmitting boundaries are important for modeling the wave propagation in the finite element analysis of dynamic foundation problems. In this study, viscoelastic boundaries for multiple seismic waves or excitations sources were derived for two-dimensional and three-dimensional conditions in the time domain, which were proved to be solid by finite element models. Then, the method for equivalent forces’ input of seismic waves was also described when the proposed artificial boundaries were applied. Comparisons between numerical calculations and analytical results validate this seismic excitation input method. The seismic response of subway station under different seismic loads input methods indicates that asymmetric input seismic loads would cause different deformations from the symmetric input seismic loads, and whether it would increase or decrease the seismic response depends on the parameters of the specific structure and surrounding soil.
1. Introduction
The 1995 Kobe earthquake caused a major collapse of the Daikai subway station in Kobe, which represents the first modern underground structure failure during a seismic event. The seismic design and analysis methods of underground structures have been rapidly developed since then [1–8]. Hashash et al. [2] described approaches used by engineers in quantifying the seismic effects on an underground structure. Both analytical methods and numerical calculations approaches were reviewed in their literature. Compared with the analytical methods, the numerical analysis has advantages in handling the problems of soil-structure interaction, inelastic behavior of the structure, and surrounding soil and the complex geometric shapes of target structures. Hence, numerical analysis methods, such as finite element methods, are widely used in the seismic analysis of underground structures.
To avoid the undesirable reflection of waves at the artificial boundaries of the finite element model, absorbing boundaries or transmitting boundaries should be applied. Problems involving inelastic behavior of the soil in the near field are most readily solved in the time domain. Wolf [9] has reviewed the commonly used time-domain transmitting boundaries. One of the most popular transmitting boundaries is the viscous boundary proposed by Lysmer and Kuhelmeyer [10], which replaces the far field with viscous damping. However, the viscous boundary would induce an undesirable drift displacement in the whole model owing to the lack of stiffness. Deeks and Randouph [11] derived a viscoelastic boundary for two-dimensional conditions based on an approximation of the form of the outward travelling waves, which eliminate the undesirable drift displacements. Liu [12] and Du [13] proposed a viscoelastic boundary for the three-dimensional conditions using the theory of elastic wave motions. Li and Song [14] proposed a general viscous-spring transmitting boundary for dynamic analysis of saturated poroelastic media based on u-U formulation of Biot equation in cylindrical coordinate. However, these viscelastic boundary conditions are derived in the media where only one excitation source exists. In fact, multiple waves or excitations sources also exist.
In this study, a viscoelastic boundary for multiple waves or excitation sources is presented and validated. The main work of this study is as follows. First, a two-dimensional viscoelastic boundary for multiple waves or excitation sources is derived based on an approximation of the form of the outward traveling waves. The proposed visc-elastic boundary is validated by a finite element model. Next, the proposed two-dimensional viscoelastic boundary is promoted for three-dimensional condition. Then, input equivalent forces of seismic waves are described when the viscoelastic boundary for multiple waves or excitation sources is used. Finally, the seismic response of a subway station using different seismic loads input methods is discussed.
2. Derivation and Validation of the Two-Dimensional Viscoelastic Boundary Conditions
In this section, a two-dimensional viscoelastic boundary for the multiple input waves or excitation sources would be developed based on an approximation of the form of the outward traveling waves.
2.1. Shear Boundary Equations
According to the elastic wave propagation theory, the two-dimensional shear waves in the elastic media followwhere represents the shear wave velocity in the media, correlated to the shear modulus and the density of the media as follows:
According to Whitham [15], the solution of (1) closely approximates to the following expression:
Assuming two independent outgoing waves and , the displacement perpendicular to the wave front of boundary point X is
The velocity of point X is as follows:
The corresponding shear strain and shear stress are expressed as in
Assume that the shear stress at point X is related to the velocity and displacement at point X as
Combining (4), (5), (7), and (8) yields the following:
An exact solution of (9) cannot be found. However, (9) could have a solution if the following were solid:
Equation (10) would be solid only when the shear stiffness parameter follows and , which yields r1 = r2. This is unavailable for most conditions. As in (11), damping parameter of the two-dimensional viscoelastic boundary conditions for multiple seismic excitation sources is same as those of viscous boundary conditions and viscoelastic boundary conditions of Deeks and Randoph [11]. The shear stiffness parameter K, however, has no accurate numerical solutions. To obtain the shear stiffness parameter for two-dimensional viscoelastic boundary conditions under multiple seismic excitation sources input, three kinds of shear stiffness parameters are proposed.
Assuming r1 is shorter than r2, if shear stiffness parameters are mainly affected by the source at distance of r1, the shear stiffness parameter K1 follows
If shear stiffness parameters are mainly affected by the source at distance of r2, the shear stiffness parameter K2 is given by
If shear stiffness parameters are affected by the both sources, the shear stiffness parameter K3 is expressed as
To employ the proposed viscoelastic boundary equations in the finite element analysis model, discrete springs and dashpots should be implemented on the corresponding boundary nodes, which would make the shear stresses and displacements approximate to the real ones. Since (8) is in the stress condition, the proposed shear stiffness and damping should be multiplied by the effective area of the meshes.
2.2. Normal Boundary Equations
As in Section 2.1, the normal stiffness parameters for two-dimensional viscoelastic boundary conditions are assumed as follows:
As in the viscous boundary, the damping parameter follows
2.3. Numerical Verification of the Two-Dimensional Viscoelastic Boundary Conditions
A homogeneous elastic soli model is constructed to validate the viscoelastic boundary conditions. Shown in Figure 1, the model (M1) is 200 meters wide and 100 meters high. The shear wave velocity is 224 m/s, density is 2000 kg/m3, and Poisson’s ratio is 0.25. A triangular impulsive excitation is imposed at sites B and C on the top of the model. The input triangular impulsive excitation is displayed in Figure 2. The other three boundaries of the model are set as different kinds of boundary conditions, including free boundary condition, solid boundary, viscous boundary, and the three proposed viscoelastic boundary conditions. Horizontal and vertical displacements at site A are recorded when different boundary conditions are implemented during the dynamic calculation. A model (M2) with larger size is constructed as a reference, shown in Figure 1. Since the reflected input triangular wave takes more than 3 seconds to site A in M2 and the displacements at site A in the first 2 seconds are considered, the reflected waves in M2 has no effects on the response of site A in this case. The displacements at site A in M2 are used as a reference to the responses of M1 under different boundary conditions.
(a)
(b)
Figure 3 shows the displacements of A subjected to the triangular excitation under different boundary conditions. In Figure 3, FB, RB, and VB denote the free surface boundary condition, rigid boundary condition, and viscous boundary condition, respectively. EV1, EV2, and EV3 represent the stiffness of viscoelastic boundary controlled by the near excitation source, the further excitation source, and both sources, respectively.
(a)
(b)
(c)
(d)
In Figures 3(a) and 3(c), horizontal and vertical displacements of site A under the free boundary, compared with the results of M2, both have an undesirable drift deformation. This is because the free surface boundary has no stiffness to recover from elastic deformation and no damping to absorb the energy of transmitting waves. Since the rigid boundary has infinitely large stiffness to recover from elastic deformation but no damping to absorb the transmitting waves, oscillation appears in M1 when rigid boundary is employed. Displacements at site A under viscous boundary approximate to the results of M2 within the first 0.5 second but show some slight drift in the following 1.5 seconds due to the lack of stiffness.
In Figures 3(b) and 3(d), also compared with M2, displacements at site A under the proposed three viscoelastic boundary conditions (EV1, EV2, and EV3) are more accurate than those under viscous boundary conditions, which eliminate the drift displacement in viscous boundary conditions. In the first 0.6 second, displacements at site A under the proposed three viscoelastic boundary conditions (EV1, EV2, and EV3) are the same. In the last 1.4 seconds, displacement waveform of EV3 is between those of EV1 and EV2.
From the above calculated results and analysis and considering that the distances between the excitation sources and boundaries could be distinctly different, it is suggested to estimate the shear stiffness parameter K by the both sources.
3. The 3D Viscoelastic Boundary Conditions
Similar to the proposed 2D viscoelastic boundary equations in Section 2, the proposed two-dimensional viscoelastic boundary is extended to the three-dimensional condition. The three-dimensional viscoelastic boundary equations proposed by Du and Zhao [13] is adopted as a reference. The shear stiffness parameter K for multiple excitation sources could be obtained bywhere n is the number of excitations. The normal stiffness parameter K for multiple excitation sources is expressed as
The damping parameters are same as the two-dimensional conditions. To verify the proposed 3D boundary condition for multiple excitation sources, a three-dimensional homogeneous elastic model is constructed with a size of , shown in Figure 4. The shear wave velocity, density, and Poisson’s ratio are 1, 1, and 0.3, respectively. The triangular impulsive excitation in Figure 2 is imposed along three lines (EE, DD, and FF) on the top of free surface. In lines DD, EE, and FF, the excitations are along x direction, y direction, and z direction, respectively. The vertical distances between center point O on the free surface and the three lines are all 0.2. Viscous boundary conditions and the suggested EV3 boundary equations are employed in M3, respectively. As a reference, a model 4 (M4) with a size of is constructed. The displacements of site O are calculated under different boundary conditions when the free surface is excited by the triangular excitation. The comparison is displayed in Figure 5.
(a)
(b)
(c)
The calculated results in three-dimensional model are similar to those in two-dimensional model. Using the results of M4FB as a reference, responses of M3EV3 are more accurate than those of M3VB. The proposed viscoelastic boundary could allow the model to recover from elastic deformation and absorb the transmitting waves. Thus, the proposed viscoelastic boundary conditions could be employed for the dynamic calculation when multiple excitation sources are imposed.
4. Input Equivalent Forces of Seismic Waves
The appropriate method for input of seismic waves is important in soil-structure interaction dynamic analysis. For different artificial boundary conditions, different methods for seismic waves input should be adopted. Joyner and Chen [16] used equivalent loads for seismic waves input in one-dimensional model when viscous boundary is employed. In this section, the equivalent loads for viscoelastic boundary under multiple excitation sources are presented.
4.1. The Method for Input of Seismic Waves
To achieve accurate input of seismic waves, the input equivalent loads should make the displacements and stress conditions of artificial boundary same as the real ones. Assuming is the original vertical shear waves on artificial boundaries and is induced displacements by input shear waves , they should followwhere is the original stress caused by
To employ the proposed viscoelastic boundary conditions, discrete springs and dashpots are implemented. is the input stress at artificial boundary site B, and is the internal stress between the elements and artificial boundary. Then, the stress of site B is calculated by
Internal stress could be obtained by
The input loads on boundary could be estimated from (22) to (25)
For the two-dimensional viscoelastic boundary under multiple excitation sources, (26) could be rewritten as follows by (11) and (14):
in (27) is the input shear stress on boundary nodes, which should be multiplied by the area of corresponding mesh during dynamic calculation. For input forces on boundary nodes in normal direction, primary wave speed should replace the shear wave speed.
For the three-dimensional viscoelastic boundary under multiple excitation sources, (26) and corresponding stiffness parameter and damping parameter could be used to calculate the input loads. Liu et al. [17] discussed the sign and direction for the three terms in (26) on the five artificial boundaries, which are adopted as a reference in this study and summarized in Table 1.
4.2. Validation of Input Equivalent Forces of Seismic Waves
To validate the proposed method for input equivalent forces of seismic waves, a homogeneous elastic model of 20 m × 20 m × 20 m is constructed. Except for the free surface, the other five artificial faces are modeled as proposed viscoelastic boundaries, shown in Figure 6. The density of soil media is 2 × 103 kg/m3, Young’s modulus is 3 MPa, shear wave velocity is 25 m/s, and Poisson’s ratio is 0.25.
Figure 7 displays the displacement history of the vertical incident shear wave. Two different methods for equivalent forces input are compared. In model 5 (5M), the equivalent forces by (26) and Table 1 are imposed on the five viscoelastic boundaries. In model 6 (1M), equivalent forces by (26) are only imposed on the bottom of the model. The displacements O at the center of free surface are calculated and compared with the analytical solution.
Figure 8(a) shows that displacements of O, calculated by viscoelastic boundary under multiple excitation sources and corresponding equivalent seismic loads, highly satisfy the analytical ones. In Figure 8(b), displacement responses at O are similar to the analytical one. But it has a peak value of 0.09 m, which is only a half of the analytical one. The discrepancy is caused by inappropriate input of seismic loads. During the incident wave propagation, boundary nodes on the four vertical boundaries would induce spring forces and damping forces due to the displacement and velocity of boundary nodes. In Figure 8(b), equivalent forces by (26) are only imposed on the bottom of the model. The stress conditions on the four vertical boundaries are different from real ones, which causes the discrepancy between the results of 1M and the analytical ones.
(a)
(b)
Comparisons between numerical and analytical results validate the proposed viscoelastic boundary under multiple excitation sources and corresponding equivalent seismic loads input method. If equivalent forces are only imposed on the bottom of the model, it would underestimate the numerical results.
5. Seismic Response of Subway Station under Different Seismic Loads Input Methods
For dynamic analysis of underground structures, appropriate equivalent seismic loads input methods should be used for different artificial boundary conditions. In Sections 2–4, the accuracy of proposed viscoelastic boundary and corresponding equivalent seismic loads input method have been validated. In real seismic events, the ground motions would vary at different sites when the structure has a large size. The difference of ground motions could induce the difference of stress condition and motion on different parts of the structure, especially for the structure with a large size. In the following section, the seismic response of a subway station subjected to different seismic waves will be analyzed. A group of different seismic waves are used as incident loads.
5.1. Analysis Model and Parameters
A 2-storey and 4-span subway station is selected as the target structure, with a cross section displayed in Figure 9. The whole model is 252 meters wide, 68 meters high, and 200 meters long. The depth of the roof is 3 meters. The roof and bottom of the station are both 0.8 meter thick. The floor 2 is 0.4 meter thick. The two side walls are both 0.8 meter thick. Center columns have a radius of 0.5 meter. The heights of storey 1 and storey 2 are 6 meters and 7 meters, respectively. Distance between two adjacent spans is 7 meters in lateral direction and 10 meters in longitudinal direction. The proposed viscoelastic boundary is employed.
Both subway station and surrounding soil are assumed to be elastic in the following analysis. Velocities of shear wave and primary wave speed in the surrounding soil are 200 m/s and 350 m/s, respectively. The detailed parameters of subway station and surrounding soil are summarized in Table 2. No slip is allowed between the station and surrounding soil during the seismic analysis.
The EW components of MXN wave in Wenchuan earthquake are selected as the input seismic wave in lateral direction, the peak acceleration of which is amplified to 5 m/s2. Seismic response analyses of the subway station under three different seismic loads input conditions are conducted. Table 3 summarizes the imposed waves on the boundaries for the three cases. In case 1, equivalent seismic forces of MXN wave are only imposed on the bottom of model. In case 2, bottom and the four vertical boundaries of model are all subjected to the corresponding equivalent seismic forces of MXN wave. In case 3, the corresponding equivalent seismic forces of MXN wave are imposed to the bottom and left side boundary of the model; corresponding equivalent seismic forces of SFB wave, LDD wave, and JYC wave are imposed to the right side, front, and back boundary of the model, respectively. In these three cases, the peak acceleration of incident waves is all amplified to 5 m/s2, as displayed in Figure 10.
(a) JYC wave
(b) LDD wave
(c) MXN wave
(d) SFB wave
5.2. Dynamic Calculation and Analyses
Story drift displacement histories of the floor 1 are compared among different cases, as shown in Figure 11. The waveforms of the three cases are similar, but the amplitude differs from each other. In Figure 11(a), the amplitudes of case 1 are smaller than those of case 2, which is caused by the difference of seismic forces input methods. In Figure 11(b), the amplitudes of case 3 are larger than those of case 2, which is induced by the difference of incident waves. Although the peak accelerations of four incident waves are all amplified to 5 m/s in case 3, the four waves have different spectral components in frequency domain.
(a) case 1 and case 2
(b) case 2 and case 3
Figure 12 displays the maximum story drift displacements at different parts of the station for three cases. The results show that maximum story drift displacements in case 1 are only 67% of those in case 2. As explained in Section 4.2, this discrepancy is caused by the difference of seismic forces input methods. If equivalent seismic forces are only imposed on the bottom of model, maximum story drift displacements of the station would be underestimated. In case 1 and case 2, the maximum story drift displacements of left columns approximate to those of right columns, which are slightly larger than those of center columns. In case 3, the maximum story drift displacements of right columns are distinctly larger than those of left columns. In case 3, the distribution of maximum story drift displacements changes, due to the asymmetric input seismic forces, while the input forces are both symmetric in case 1 and case 2, which induce symmetric maximum story drift deformation.
Variable equivalent seismic forces on vertical boundaries would not always induce a larger story drift displacement than that induced by uniform equivalent seismic forces on vertical boundaries. The seismic analysis in case 3 could only indicate that asymmetric seismic forces input would cause asymmetric deformation. Compared to the symmetric seismic forces input, whether it would increase or decrease the seismic response depends on the parameters of the specific structure and surrounding soil.
6. Conclusions
In this study, the viscoelastic boundaries for multiple seismic waves or excitations sources have been derived and validated for two-dimensional and three-dimensional conditions in the time domain. The corresponding method for equivalent seismic loads input has also been presented.
First, a two-dimensional viscoelastic boundary for multiple waves or excitation sources was derived based on an approximation of the form of the outward traveling waves. The proposed two-dimensional viscoelastic boundary was promoted to the three-dimensional condition. The proposed viscoelastic boundary was validated using a finite element model by the numerical calculation.
Then, the method for input equivalent forces of seismic waves was described when the viscoelastic boundary for multiple waves or excitation sources was used.
Finally, the seismic responses of a subway station using different seismic loads input methods were discussed. Compared with the symmetric seismic forces input, whether the asymmetric seismic forces input would increase or decrease the seismic response depends on the parameters of the specific structure and surrounding soil.
Data Availability
All the data used to support the findings of this study are included within the article. Requests for data are available from the corresponding author after this article is published.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (Grant no. 51708018) and Beijing Postdoctoral Science Foundation (Grant no. 2017ZZ063).