# Viscoelastic Boundary Conditions for Multiple Excitation Sources in the Time Domain

Authors: Chen Xia; Chengzhi Qi; Xiaozhao Li
Journal: Mathematical Problems in Engineering
 (2018)

---

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


---

## Body

## 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 follow(1)∂2u∂t2=cs2∂2u∂r2+1r∂u∂r,where cs represents the shear wave velocity in the media, correlated to the shear modulus G and the density ρ of the media as follows:(2)cs=Gρ.According to Whitham [15], the solution of (1) closely approximates to the following expression:(3)ur,t=1rfr-cst.Assuming two independent outgoing wavesu1(r1,t) and u2(r2,t), the displacement perpendicular to the wave front of boundary point X is(4)uX=u1r1,t+u2r2,t.The velocity of point X is as follows:(5)uX•=-csr1f′r1-cst-csr2f′r2-cst.The corresponding shear strain and shear stress are expressed as in(6)γX=-∂uX∂r=12r11.5fr1-cst-1r1f′r1-cst+12r21.5fr2-cst-1r2f′r2-cst,(7)τX=G12r11.5fr1-cst-1r1f′r1-cst+12r21.5fr2-cst-1r2f′r2-cst.Assume that the shear stress at point X is related to the velocity and displacement at point X as(8)τX=KuX+CuX•.Combining (4), (5), (7), and (8) yields the following:(9)G12r11.5fr1-cst-1r1f′r1-cst+12r21.5fr2-cst-1r2f′r2-cst=Kr1fr1-cst+Kr2fr2-cst-Ccsr1f′r1-cst-Ccsr2f′r2-cst.An exact solution of (9) cannot be found. However, (9) could have a solution if the following were solid:(10)G2r11.5fr1-cst+G2r21.5fr2-cst=Kr1fr1-cst+Kr2fr2-cst,(11)C=Gcs.Equation (10) would be solid only when the shear stiffness parameter K follows K=G/2r1 and K=G/2r2, which yieldsr1 =r2. This is unavailable for most conditions. As in (11), damping parameter C 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 parameterK, 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.Assumingr1 is shorter thanr2, if shear stiffness parameters are mainly affected by the source at distance ofr1, the shear stiffness parameterK1 follows(12)K1=G2r1.If shear stiffness parameters are mainly affected by the source at distance ofr2, the shear stiffness parameterK2 is given by(13)K2=G2r2.If shear stiffness parameters are affected by the both sources, the shear stiffness parameterK3 is expressed as(14)K3=G2r1r2.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 Section2.1, the normal stiffness parameters for two-dimensional viscoelastic boundary conditions are assumed as follows:(15)K1=2Gr1,(16)K2=2Gr2,(17)K3=2Gr1r2.As in the viscous boundary, the damping parameterC follows(18)C=ρcp.


### 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 Figure1, 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.Figure 1
The 2D model of rectangular soil rectangular: a represents M1 and b represents M2.
(a) (b)Figure 2
Input triangular loads history at points B and C.Figure3 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.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.
(a) (b) (c) (d)In Figures3(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 Figures3(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 parameterK by the both sources.


## 2.1. Shear Boundary Equations

According to the elastic wave propagation theory, the two-dimensional shear waves in the elastic media follow(1)∂2u∂t2=cs2∂2u∂r2+1r∂u∂r,where cs represents the shear wave velocity in the media, correlated to the shear modulus G and the density ρ of the media as follows:(2)cs=Gρ.According to Whitham [15], the solution of (1) closely approximates to the following expression:(3)ur,t=1rfr-cst.Assuming two independent outgoing wavesu1(r1,t) and u2(r2,t), the displacement perpendicular to the wave front of boundary point X is(4)uX=u1r1,t+u2r2,t.The velocity of point X is as follows:(5)uX•=-csr1f′r1-cst-csr2f′r2-cst.The corresponding shear strain and shear stress are expressed as in(6)γX=-∂uX∂r=12r11.5fr1-cst-1r1f′r1-cst+12r21.5fr2-cst-1r2f′r2-cst,(7)τX=G12r11.5fr1-cst-1r1f′r1-cst+12r21.5fr2-cst-1r2f′r2-cst.Assume that the shear stress at point X is related to the velocity and displacement at point X as(8)τX=KuX+CuX•.Combining (4), (5), (7), and (8) yields the following:(9)G12r11.5fr1-cst-1r1f′r1-cst+12r21.5fr2-cst-1r2f′r2-cst=Kr1fr1-cst+Kr2fr2-cst-Ccsr1f′r1-cst-Ccsr2f′r2-cst.An exact solution of (9) cannot be found. However, (9) could have a solution if the following were solid:(10)G2r11.5fr1-cst+G2r21.5fr2-cst=Kr1fr1-cst+Kr2fr2-cst,(11)C=Gcs.Equation (10) would be solid only when the shear stiffness parameter K follows K=G/2r1 and K=G/2r2, which yieldsr1 =r2. This is unavailable for most conditions. As in (11), damping parameter C 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 parameterK, 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.Assumingr1 is shorter thanr2, if shear stiffness parameters are mainly affected by the source at distance ofr1, the shear stiffness parameterK1 follows(12)K1=G2r1.If shear stiffness parameters are mainly affected by the source at distance ofr2, the shear stiffness parameterK2 is given by(13)K2=G2r2.If shear stiffness parameters are affected by the both sources, the shear stiffness parameterK3 is expressed as(14)K3=G2r1r2.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 Section2.1, the normal stiffness parameters for two-dimensional viscoelastic boundary conditions are assumed as follows:(15)K1=2Gr1,(16)K2=2Gr2,(17)K3=2Gr1r2.As in the viscous boundary, the damping parameterC follows(18)C=ρcp.


## 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 Figure1, 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.Figure 1
The 2D model of rectangular soil rectangular: a represents M1 and b represents M2.
(a) (b)Figure 2
Input triangular loads history at points B and C.Figure3 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.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.
(a) (b) (c) (d)In Figures3(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 Figures3(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 parameterK by the both sources.


## 3. The 3D Viscoelastic Boundary Conditions

Similar to the proposed 2D viscoelastic boundary equations in Section2, 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 parameterK for multiple excitation sources could be obtained by(19)K=2Gr1•r2•••rnn,wheren is the number of excitations. The normal stiffness parameterK for multiple excitation sources is expressed as(20)K=4Gr1•r2•••rnn.The damping parametersC 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.Figure 4
Locations of input seismic loads and observation site on the free surface of the 3D model.Figure 5
a represents the displacements of O in the x direction for different model size and artificial boundary conditions, b represents the displacements in the y direction, and c represents displacements of O in the z direction.
(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. Assumingu0(xB,yB,t) is the original vertical shear waves on artificial boundaries and u(xB,yB,t) is induced displacements by input shear waves τ(xB,yB,t), they should follow(21)uxB,yB,t=u0xB,yB,t,(22)τxB,yB,t=τ0xB,yB,t,where τ0(xB,yB,t) is the original stress caused by u0(xB,yB,t)(23)τ0xB,yB,t=ρcsu•0xB,yB,t.To employ the proposed viscoelastic boundary conditions, discrete springs and dashpots are implemented.FB(t) is the input stress at artificial boundary site B, and fB(t) is the internal stress between the elements and artificial boundary. Then, the stress of site B is calculated by(24)τxB,yB,t=FBt-fBt.Internal stressfB(t) could be obtained by(25)fBt=Ku0xB,yB,t+Cu•0xB,yB,t.The input loads on boundary could be estimated from (22) to (25)(26)FBt=Ku0xB,yB,t+Cu•0xB,yB,t+ρcsu•0xB,yB,t.For the two-dimensional viscoelastic boundary under multiple excitation sources, (26) could be rewritten as follows by (11) and (14):(27)FBt=G2r1r2u0xB,yB,t+2ρcsu•0xB,yB,t.F 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 cp should replace the shear wave speed.For the three-dimensional viscoelastic boundary under multiple excitation sources, (26) and corresponding stiffness parameter K and damping parameter C 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.Table 1
Parameters of the subway station and surrounding soil.
Boundary Left Right Front Back Bottom Input in the x direction Stress term z+ z- x+ x+ x+ Velocity and displacement term x+ x+ - - x+ Input in the y direction Stress term z+ z- - - y+ Velocity and displacement term y+ y+ - - y+ Input in the z direction Stress term x- x+ y- y+ z- Velocity and displacement term z- z- z- z- z- 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.


### 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 Figure6. 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 6
Three-dimensional model of the soil and viscoelastic artificial boundary condition.Figure7 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 7
Input sinusoidal displacement waveform.Figure8(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.Figure 8
The displacement waveforms of observation point O under different conditions.
(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.


## 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. Assumingu0(xB,yB,t) is the original vertical shear waves on artificial boundaries and u(xB,yB,t) is induced displacements by input shear waves τ(xB,yB,t), they should follow(21)uxB,yB,t=u0xB,yB,t,(22)τxB,yB,t=τ0xB,yB,t,where τ0(xB,yB,t) is the original stress caused by u0(xB,yB,t)(23)τ0xB,yB,t=ρcsu•0xB,yB,t.To employ the proposed viscoelastic boundary conditions, discrete springs and dashpots are implemented.FB(t) is the input stress at artificial boundary site B, and fB(t) is the internal stress between the elements and artificial boundary. Then, the stress of site B is calculated by(24)τxB,yB,t=FBt-fBt.Internal stressfB(t) could be obtained by(25)fBt=Ku0xB,yB,t+Cu•0xB,yB,t.The input loads on boundary could be estimated from (22) to (25)(26)FBt=Ku0xB,yB,t+Cu•0xB,yB,t+ρcsu•0xB,yB,t.For the two-dimensional viscoelastic boundary under multiple excitation sources, (26) could be rewritten as follows by (11) and (14):(27)FBt=G2r1r2u0xB,yB,t+2ρcsu•0xB,yB,t.F 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 cp should replace the shear wave speed.For the three-dimensional viscoelastic boundary under multiple excitation sources, (26) and corresponding stiffness parameter K and damping parameter C 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.Table 1
Parameters of the subway station and surrounding soil.
Boundary Left Right Front Back Bottom Input in the x direction Stress term z+ z- x+ x+ x+ Velocity and displacement term x+ x+ - - x+ Input in the y direction Stress term z+ z- - - y+ Velocity and displacement term y+ y+ - - y+ Input in the z direction Stress term x- x+ y- y+ z- Velocity and displacement term z- z- z- z- z- 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.


## 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 Figure6. 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 6
Three-dimensional model of the soil and viscoelastic artificial boundary condition.Figure7 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 7
Input sinusoidal displacement waveform.Figure8(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.Figure 8
The displacement waveforms of observation point O under different conditions.
(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 Figure9. 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.Figure 9
The cross section of the subway station.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 Table2. No slip is allowed between the station and surrounding soil during the seismic analysis.Table 2
Parameters of the subway station and surrounding soil.
Density(kg/m3) Young’s modulus (MPa) Poisson’s ratio Station 2.5 × 103 2.8 × 104 0.20 Soil 2.0 × 103 3.0 × 102 0.25The 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.Table 3
The imposed waves on the boundaries for the three cases.
Case Bottom Left Right Front Back 1 MXN - - - - 2 MXN MXN MXN MXN MXN 3 MXN MXN SFB LDD JYCFigure 10
The selected acceleration waveforms of Wenchuan earthquake.
(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 Figure11. 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.Figure 11
Story drift displacements of left floor 1 columns in the subway station under different conditions.
(a)
case 1 and case 2 (b)
case 2 and case 3Figure12 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.Figure 12
Maximum story drift angels of columns in the subway station for different cases.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.


## 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 Figure9. 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.Figure 9
The cross section of the subway station.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 Table2. No slip is allowed between the station and surrounding soil during the seismic analysis.Table 2
Parameters of the subway station and surrounding soil.
Density(kg/m3) Young’s modulus (MPa) Poisson’s ratio Station 2.5 × 103 2.8 × 104 0.20 Soil 2.0 × 103 3.0 × 102 0.25The 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.Table 3
The imposed waves on the boundaries for the three cases.
Case Bottom Left Right Front Back 1 MXN - - - - 2 MXN MXN MXN MXN MXN 3 MXN MXN SFB LDD JYCFigure 10
The selected acceleration waveforms of Wenchuan earthquake.
(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 Figure11. 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.Figure 11
Story drift displacements of left floor 1 columns in the subway station under different conditions.
(a)
case 1 and case 2 (b)
case 2 and case 3Figure12 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.Figure 12
Maximum story drift angels of columns in the subway station for different cases.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.



---

*Source: 7982342-2018-10-10.xml*