Viscoelastic Boundary Conditions for Multiple Excitation Sources in the Time Domain

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 describedwhen 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.


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][2][3][4][5][6][7][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 twodimensional viscoelastic boundary is promoted for threedimensional 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.

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.

Shear Boundary Equations.
According to the elastic wave propagation theory, the two-dimensional shear waves in the elastic media follow where  s 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  1 ( 1 , ) and  2 ( 2 , ), 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  = /2 1 and  = /2 2 , which yields r 1 = r 2 .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 r 1 is shorter than r 2 , if shear stiffness parameters are mainly affected by the source at distance of r 1 , the shear stiffness parameter K 1 follows If shear stiffness parameters are mainly affected by the source at distance of r 2 , the shear stiffness parameter K 2 is given by If shear stiffness parameters are affected by the both sources, the shear stiffness parameter K 3 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.

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  =  푝 .
(18) 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.

Numerical Verification of the Two-Dimensional Viscoelastic Boundary
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.

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 by where 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 1 × 1 × 0.5, 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 4 × 4 × 2 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.
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.

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 onedimensional model when viscous boundary is employed.In this section, the equivalent loads for viscoelastic boundary under multiple excitation sources are presented.

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  0 ( 퐵 ,  퐵 , ) is the original vertical shear waves on artificial boundaries and ( 퐵 ,  퐵 , ) is induced displacements by input shear waves ( 퐵 ,  퐵 , ), they should follow where  0 ( 퐵 ,  퐵 , ) is the original stress caused by  0 ( 퐵 ,  퐵 , ) To employ the proposed viscoelastic boundary conditions, discrete springs and dashpots are implemented. B (t) is the input stress at artificial boundary site B, and  B (t) is the internal stress between the elements and artificial boundary.Then, the stress of site B is calculated by Internal stress  B (t) 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): B (t) 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.

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 × 10 3 kg/m 3 , 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  The letter "x+" in the table means that the corresponding term is in the positive x direction and other letters follow the same rule.The minus signs "-" means that there are no corresponding terms for these boundaries.Stress, velocity, and displacement terms are the third, second, and first term in (26), respectively.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.

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.and 350 m/s, respectively.The detailed parameters of subway station and surrounding soil are summarized in Table 2.

Analysis Model and
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/s 2 .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 Table 3: The imposed waves on the boundaries for the three cases.

Case
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/s 2 , as displayed in Figure 10. 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.

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.
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/m 3 , 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 Figure1.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.

Figure 1 :Figure 2 :
Figure 1: The 2D model of rectangular soil rectangular: a represents M1 and b represents M2.

Figure 3 :
Figure 3: a and b represent the horizontal displacements at site A for different model size and artificial boundary conditions, while c and d represent the vertical displacements.

Figure 4 :Figure 5 :
Figure 4: Locations of input seismic loads and observation site on the free surface of the 3D model.

Figure 6 :
Figure 6: Three-dimensional model of the soil and viscoelastic artificial boundary condition.

Figure 8 :Figure 9 :
Figure 8: The displacement waveforms of observation point O under different conditions.

5. 2 .Figure 12 :
Figure 12: Maximum story drift angels of columns in the subway station for different cases.

Table 1 :
Parameters of the subway station and surrounding soil.