# Numerical Analysis on Seismic Responses of a Large Metro Station in Sandy Area

Authors: Yu Liu; Xiao Liu
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

A new assessment method is proposed considering shortcomings of the current investigations for the assessment of underground structures considering seismic load. First, the dynamic elastic modulus and damping ratio of sand are studied based on the dynamic characteristic experiment of the soil. Subsequently, a series of design acceleration time history curves are obtained using the proposed methods. Finally, the seismic response of underground structures is evaluated by a three-dimensional numerical simulation. The results indicate that the dynamic elastic modulus of sand decreases with the increase of dynamic elastic strain and increases with the increase of confining pressure. The influence of confining pressure on the damping ratio of sand is slight, but it can still be seen that the damping ratio decreases with the increase of confining pressure, especially when the dynamic elastic strain is relatively low. The interlayer displacement between the top and middle plate is 2.1 mm and the interlayer displacement angle is 1/2650. The interlayer displacement between the middle and bottom plate is 1.7 mm and the interlayer displacement angle is 1/4117. Under the action of fortification earthquake, the interlayer displacement angles of Taiyuan Street station are less than 1/550, indicating that the structure meets the seismic requirements.


---

## Body

## 1. Introduction

With the strategy of a country with a strong transportation network, the construction of an urban subway has entered a stage of rapid development in China. Although the complex structure of the subway station has rich functions and can promote regional development, its seismic response is also difficult to obtain accurately because of its complex structure under earthquakes [1, 2]. With the emergence of cases of serious damage to underground structures caused by earthquakes, more and more scholars began to do a lot of research work in this field [3, 4].In the theoretical analysis, Newmark [5] and Kuesel [6] proposed a simplified algorithm for free field strain under the condition of the simple harmonic incident at any angle in isotropic elastic homogenization. Luco and Barros [7] studied the seismic response of circular underground structures in semi-infinite elastic space under SH wave by using the method of combining finite elements with a wave function and the indirect boundary integral method based on two-dimensional Green’s function. Davis et al. [8] derived the analytical solution of P-wave and SV-wave scattering by underground structures in semi-infinite space through the large arc assumption method. Pang et al. [9] proposed a novel model to analyze the seismic response and reliability level of a subway station, and the presented model can generate completely nonstationary ground motions. The characteristics of earthquakes have notable influences on the dynamic structural behavior of underground structures. Qiu et al. [10] proposed a modified simplified analysis method to evaluate seismic responses of a large subway station in a complex area. Based on the presented model, the seismic response of the support structure was divided into two categories such as the kinematic interaction effect and the inertial interaction effect.In the physical model experiments, Tomari and Towhata [11] used a shaking table test to explore the seismic response of flexible section structure under the condition of a liquefiable site and discussed the influence of natural vibration frequency and backfill dilatancy on structural characteristics. The shaking tests on the seismic response characteristics of the underground structure in the liquefiable site are investigated by An et al. [12]. Wu et al. [13] studied the dynamic response and spectrum characteristics of a tunnel using shaking table tests. A centrifuge shaker model test of a subway station on a liquefiable foundation was carried out, which studied the process in which damage occurs during an earthquake. Several dynamic centrifuge tests were conducted to investigate the seismic structural response and failure model of the underground structure system. The vertical load notably increases the axial forces on the frame columns to reduce the horizontal deformation capacity [14]. Miao et al. [15] investigated soil-structure interactions on the seismic response of metro stations using a series of experiments.In the aspect of numerical simulation, Chen [16] et al. conducted numerical research on the dynamic response of subway stations under seismic load in shallow soft soil sites. Taking the interval tunnel of Beijing Metro Line 7 as an example, Li et al. [17] established the numerical model of the cross tunnel by using FLAC3D software and studied the seismic response characteristics of a closed pasted cross tunnel under strong earthquakes. Lu and Huang [18] established the three-dimensional refined model of the upper and lower intersecting metro tunnel using MIDAS software, and the dynamic response law of the subway cross tunnel is analyzed. Ma et al. [19] established the finite element model of tunnel soil and the dynamic response of the silty sand layer is analyzed. Liu et al. [20] proposed a seismic control technology for shallow buried underground structures. The finite element models are established and evaluated for the support system. The results indicate that the technology can effectively reduce the horizontal displacement of the central support system and provide a reference for the ability of earthquake prevention. Du et al. [21, 22] proved that the axial compression ratio of the central column increased notably under the combined action of horizontal-vertical loads, which reduced deformation capacity in the horizontal direction.In summary, the dynamic response of a large metro station has been investigated using theoretical analyses, physical experiments, and numerical simulations. Due to the rapid development of computer technology, numerical software is widely used in geotechnical engineering. However, it is little for investigation of the seismic responses using numerical software considering the local seismic load. Thereby, considering shortcomings of the aforementioned investigations, the structural deformation of subway station under seismic load is further evaluated and several key contents are as follows: (1) The dynamic elastic modulus and damping ratio of sand are studied based on the dynamic characteristic experiment of soil. (2) A series of design acceleration time history curves are obtained using the proposed methods. (3) The seismic response of underground structures is evaluated by a three-dimensional numerical simulation.


## 2. Engineering Survey

Taiyuan Street Station of line 4 is located in a prosperous commercial area, mainly surrounded by shopping malls and business buildings. The passenger flow is relatively concentrated and the risk source is large, as shown in Figure1. The auxiliary structure of Taiyuan Street station includes 3 entrances and exits, 2 air ducts, 1 transfer channel, and 2 emergency exits. The station is a three-span island platform station on the second floor underground, with a length of 231.7 m, a standard section width of 25.3 m, and a platform width of 14 m. The buried depth of the bottom plate of the station is approximately 27 m and the covering soil of the top plate is 10.57 m. The main structure of the station is constructed by the PBA method of concealed excavation, as shown in Figure 2. The main measures of the PBA method are as follows: bored cast-in-place piles(φ1000@1200) are mainly used for side piles, a steel support is set, and the dewatering scheme outside the pit is adopted. The proposed site is mainly composed of miscellaneous fill, a cohesive soil layer, a silty soil layer, and a sandy soil layer (Table 1).Figure 1
General plan.Figure 2
Cross section.


## 3. Experimental Study on Dynamic Characteristics of the Soil

Considering that the dynamic shear modulus and damping ratio of soil are important parameters of soil dynamic characteristics, they are indispensable in the seismic safety evaluation of the engineering site and the seismic response analysis of the soil layer. The rationality of parameter selection will directly affect the safety and economy of the engineering building structure. The variation law of dynamic elastic modulus and damping ratio of sand is studied in this section..Table 1
Soil parameters.
SoilDepth (m)Weight (kN/m3)Dynamic elastic modulus (MPa)Dynamic Poisson ratioCohesion (kPa)Internal friction angle (°)Backfill soil2.01793.70.4106Fine sand2.719.5279.622326Medium coarse sand2.219.8611.530032Gravel sand8.720660.932034Silty clay2.119283.3201814Boulder clay720.01036.8——


### 3.1. Experiment Preparation

The dynamic characteristic test of soil was completed with the assistance of the laboratory of Northeast University. The indoor test was carried out according to the typical sandy soil in Shenyang. Silt is the soil sample of the bidirectional dynamic triaxial test system produced by GDS company, as shown in Figure3. Remolded soil samples with a diameter of 39.1 mm and a height of 80 mm were used in the test. The remolded sample is prepared by the multilayer wet tamping method, which is carried out in five layers. The dry density of the silty sand sample is 1.55 g/cm3. The consolidation stress of the silt sample is 50 kPa, 100 kPa, and 150 kPa respectively. The weight of each layer of the soil sample is determined according to the dry density and predesigned water content of the soil sample and compacted to the corresponding height. The contact surface of each layer is scratched to ensure a good upper and lower contact.Figure 3
Dynamic triaxial test of sandy soil.


### 3.2. Result Analysis

#### 3.2.1. Dynamic Elastic Modulus

The relationship curve between dynamic elastic modulus and dynamic elastic strain of the sample under different confining pressures is shown in Figure4. It can be seen from the figure that with the gradual increase of dynamic elastic strain, the dynamic elastic modulus decreases and the stiffness softening phenomenon occurs. At the beginning of the cycle, the curve is steep and the stiffness softening rate is fast. Then, with the increase of dynamic elastic strain, the curve tends to be flat and the stiffness softening rate decreases. For the same strain level, when the confining pressure increases, the dynamic elastic modulus increases. The dynamic elastic modulus tends to increase with the increase of confining pressure because the void ratio of the sample decreases with the increase of confining pressure, the relative density increases, and the contact point of soil particles increases, which makes the stress wave propagate faster in the soil, thus increasing the dynamic elastic modulus. The dynamic elastic modulus of sand decreases with the increase of dynamic elastic strain and increases with the increase of confining pressure.Figure 4
Ed-εd relation curve of silt under different confining pressures.


#### 3.2.2. Damping Ratio

The relationship curve between the damping ratio and dynamic elastic strain under different confining pressures is shown in Figure5. It can be seen from the figure that the damping ratio of sand increases with the increase of dynamic elastic strain. In a small strain, the damping ratio increases rapidly with the increase of dynamic elastic strain; then, the curve tends to be flat. It shows that the change of strain lags behind the change of stress in the process of vibration is limited. At the same time, the influence of confining pressure on the damping ratio of sand is not significant, but it can still be seen that the damping ratio decreases with the increase of confining pressure, especially when the dynamic elastic strain is relatively low.Figure 5
λ-εd relation curve of silt under different confining pressures.


## 3.1. Experiment Preparation

The dynamic characteristic test of soil was completed with the assistance of the laboratory of Northeast University. The indoor test was carried out according to the typical sandy soil in Shenyang. Silt is the soil sample of the bidirectional dynamic triaxial test system produced by GDS company, as shown in Figure3. Remolded soil samples with a diameter of 39.1 mm and a height of 80 mm were used in the test. The remolded sample is prepared by the multilayer wet tamping method, which is carried out in five layers. The dry density of the silty sand sample is 1.55 g/cm3. The consolidation stress of the silt sample is 50 kPa, 100 kPa, and 150 kPa respectively. The weight of each layer of the soil sample is determined according to the dry density and predesigned water content of the soil sample and compacted to the corresponding height. The contact surface of each layer is scratched to ensure a good upper and lower contact.Figure 3
Dynamic triaxial test of sandy soil.


## 3.2. Result Analysis

### 3.2.1. Dynamic Elastic Modulus

The relationship curve between dynamic elastic modulus and dynamic elastic strain of the sample under different confining pressures is shown in Figure4. It can be seen from the figure that with the gradual increase of dynamic elastic strain, the dynamic elastic modulus decreases and the stiffness softening phenomenon occurs. At the beginning of the cycle, the curve is steep and the stiffness softening rate is fast. Then, with the increase of dynamic elastic strain, the curve tends to be flat and the stiffness softening rate decreases. For the same strain level, when the confining pressure increases, the dynamic elastic modulus increases. The dynamic elastic modulus tends to increase with the increase of confining pressure because the void ratio of the sample decreases with the increase of confining pressure, the relative density increases, and the contact point of soil particles increases, which makes the stress wave propagate faster in the soil, thus increasing the dynamic elastic modulus. The dynamic elastic modulus of sand decreases with the increase of dynamic elastic strain and increases with the increase of confining pressure.Figure 4
Ed-εd relation curve of silt under different confining pressures.


### 3.2.2. Damping Ratio

The relationship curve between the damping ratio and dynamic elastic strain under different confining pressures is shown in Figure5. It can be seen from the figure that the damping ratio of sand increases with the increase of dynamic elastic strain. In a small strain, the damping ratio increases rapidly with the increase of dynamic elastic strain; then, the curve tends to be flat. It shows that the change of strain lags behind the change of stress in the process of vibration is limited. At the same time, the influence of confining pressure on the damping ratio of sand is not significant, but it can still be seen that the damping ratio decreases with the increase of confining pressure, especially when the dynamic elastic strain is relatively low.Figure 5
λ-εd relation curve of silt under different confining pressures.


## 3.2.1. Dynamic Elastic Modulus

The relationship curve between dynamic elastic modulus and dynamic elastic strain of the sample under different confining pressures is shown in Figure4. It can be seen from the figure that with the gradual increase of dynamic elastic strain, the dynamic elastic modulus decreases and the stiffness softening phenomenon occurs. At the beginning of the cycle, the curve is steep and the stiffness softening rate is fast. Then, with the increase of dynamic elastic strain, the curve tends to be flat and the stiffness softening rate decreases. For the same strain level, when the confining pressure increases, the dynamic elastic modulus increases. The dynamic elastic modulus tends to increase with the increase of confining pressure because the void ratio of the sample decreases with the increase of confining pressure, the relative density increases, and the contact point of soil particles increases, which makes the stress wave propagate faster in the soil, thus increasing the dynamic elastic modulus. The dynamic elastic modulus of sand decreases with the increase of dynamic elastic strain and increases with the increase of confining pressure.Figure 4
Ed-εd relation curve of silt under different confining pressures.


## 3.2.2. Damping Ratio

The relationship curve between the damping ratio and dynamic elastic strain under different confining pressures is shown in Figure5. It can be seen from the figure that the damping ratio of sand increases with the increase of dynamic elastic strain. In a small strain, the damping ratio increases rapidly with the increase of dynamic elastic strain; then, the curve tends to be flat. It shows that the change of strain lags behind the change of stress in the process of vibration is limited. At the same time, the influence of confining pressure on the damping ratio of sand is not significant, but it can still be seen that the damping ratio decreases with the increase of confining pressure, especially when the dynamic elastic strain is relatively low.Figure 5
λ-εd relation curve of silt under different confining pressures.


## 4. Site Design Ground Motion Parameters

### 4.1. Synthetic Ground Motion Time History of Bedrock

In the seismic response time history analysis, the selection of an appropriate ground motion acceleration time process is very important. Although the number of actual seismic records has greatly increased in the past few decades, the site conditions of the recording sites may be very different from the construction site conditions we are concerned about. We need a set of ground motion samples that meet the same statistical characteristics as the seismic response input. In this section, based on the above-ground motion characteristic parameters of the project site, the triangular series superposition simulation method is used to synthesize the bedrock ground motion time history. The calculation model is as follows:(1)x″t=ft∑k=0nCkcosωkt+φk=ftat,(2)at=∑k=0nCkcosωkt+φk.In the formula,φk is the random phase angle uniformly distributed in the 0,2π interval; Ck and ωk is the amplitude and frequency of the k-th frequency component, respectively; ft is the strength envelope function, which is a definite function of time; and at is a stationary Gaussian process.It can be seen from equation (1) that the synthetic acceleration time history x″t with the characteristics of a nonstationary random process is the product of a Gaussian process at with the characteristics of a stationary random process and an intensity envelope function ft.The coefficientCk in equation (2) can be determined by a given power spectral density function Sωk , where we obtain as follows:(3)Ck=4SωkΔω,ωk=2πkT,Δω=2πT,where t is the total holding time of the stationary random process. In order to use the acceleration response spectrum as the target spectrum of artificial acceleration time history, the following approximate relationship between the response spectrum and power spectrum can be used to replace Sωk in equation (3);(4)Sωk=−ζπωSaTω21ln−π/ωTln1−P.In the formula,SaTω is the given target acceleration response spectrum, ζ is the damping ratio, t is the duration, and P is the exceedance probability of the response.The intensity envelope functionft is related to the magnitude, epicenter distance, geological structure background, site soil conditions, and other factors and reflects the unstable characteristics of ground motion with time. Generally, the strength envelope function ft is expressed by the following piecewise function:(5)ft=t/t12,0≤t≤t1,1,t1<t≤t2,e−ct−t2,t2<t≤td,where 0,t1 is the rising section, t1,t2 is the stationary section, t2,td is the attenuation section, and C is the attenuation coefficient. Since the relationship between the response spectrum and the power spectrum represented by equation (4) is approximate, the response spectrum calculated according to the initial time history is generally only approximate to the target response spectrum. In order to improve the fitting accuracy, iterative adjustment is required. When the difference between the calculated response spectrum and the target response spectrum is less than the control accuracy (take 5%), the iteration is stopped.


### 4.2. Dynamic Response Analysis of Site Soil

#### 4.2.1. Calculation Model of Dynamic Response of the Site Soil Layer

The results of site engineering seismic conditions show that within the scope of the project site, the changes in medium characteristics and terrain are not very significant along the horizontal direction. Therefore, the influence of site conditions on seismic ground motion is considered based on a one-dimensional site model. The equivalent linearization method of one-dimensional soil shear dynamic response analysis is used in the study of the influence of the one-dimensional site model on earthquake ground motion. Its basic principle is as follows:It is assumed that the shear wave is incident vertically from the viscoelastic semi-infinite bedrock space into the horizontal layered (N-layer) nonlinear soil and propagates upward. For this calculation model, according to the wave propagation theory, the dynamic response value of the field ground medium can be calculated by using the time-frequency transformation technology, combined with the complex damping simulation of the nonlinear characteristics of soil and the equivalent linearization method.A shear harmonic is set to incident vertically upward from the calculated base and propagate in the soil layer. According to the wave theory and complex damping theory, the medium motion in each soil layer must meet the wave equation:(6)ρj∂2Ujx.t∂t2=Gjc∂2Ujx,t∂x2,where Ujx,t is the displacement value of medium reaction in the soil layer j, ρj is the density of the medium in the soil layer j, Gjc is the dynamic complex shear modulus of the medium in the soil layer j, Gjc is given by the following equation:(7)Gjc=1+2λjγjeiGjdγjeGjo,whereGjo is the maximum dynamic shear modulus of medium in soil layer, jGjdγje,λjγje is the dimensionless coefficient of equivalent dynamic shear modulus and hysteretic damping ratio of the medium in soil layer, and jγje is the maximum dynamic shear modulus of the medium in soil layer jThe medium motion between soil layers meets the conditions of displacement continuity and stress continuity.(8)Ujx,tx=Hj=Uj+1x,tx=0,τjx,tx=Hj=τj+1x,tx=0,τ1x,tx=0=0.In the formula,Hj is the layer thickness of the j-th soil layer, and it is specified that the vertically downward direction is the positive direction of the X coordinate and the coordinate origin is placed at the top of each soil layer. Through formula (6), the frequency domain value of the medium reaction quantity in the soil layer can be obtained by using the known calculated base incident wave value, and then, the time domain value of the medium reaction quantity in the soil layer can be obtained by using the Fourier transform method.Considering the nonlinear characteristics of soil, the dimensionless coefficient of equivalent dynamic shear modulus and hysteretic damping ratio of each soil layer are functions of equivalent shear strain. Therefore, in the actual calculation, the initial equivalent dynamic shear strain of the medium reaction in each soil layer is assumed and the response calculation is carried out by using the abovementioned method. Then, the maximum shear strain response of the medium at the midpoint of each soil layer is calculated. Finally, the maximum shear strain of the medium reaction at the midpoint of each soil layer multiplied by the reduction coefficient (0.65) is taken as the calculated value of the equivalent shear strain of the medium in the soil layer. We compare the equivalent shear strain used in the calculation with the equivalent dynamic shear modulus and hysteretic damping ratio corresponding to the calculated equivalent shear strain. If the relative error is less than the given allowable error (0.05), it is considered that the consideration of the nonlinear characteristics of soil meets the requirements. Otherwise, we replace the initial equivalent shear strain value with the latest calculated equivalent shear strain value and repeat the abovementioned calculation process until the relative error is less than the allowable error.


#### 4.2.2. Determination of Dynamic Parameters of the Site Soil Layer

Focal characteristics, propagation path, and local site conditions are the main factors affecting the ground motion of the site, and the influence of local site conditions is particularly prominent in a small area. Therefore, the engineering geological exploration in the proposed site and the determination of the type, stratification, and thickness of the site soil are essential work in the seismic response calculation. The analysis and determination of dynamic mechanical indexes of site soil (such as bulk density, shear wave velocity, shear modulus ratio curve, and damping ratio curve) are very important for site dynamic response analysis.In this work, the existing data are sorted out according to the requirements of soil seismic response analysis and calculation. According to the borehole data of the site, a borehole with shear wave velocity data is taken as the calculation control point for each station site and station-to-station interval. In the seismic response calculation, the medium with a shear wave velocity greater than 500 m/s is determined as the seismic input interface and the shear wave velocity of the soil below this layer is greater than 500 m/s.


### 4.3. Site Design Ground Motion Parameters

On the basis of the seismic response calculation results of the site soil layer obtained in the previous section, the design ground motion parameters of the project site in this section will be given, including the design peak ground motion acceleration and acceleration response spectrum. The design peak ground acceleration is determined by comprehensively considering the peak value of the sample and the short-period acceleration response spectrum. The design acceleration response spectrum is determined by using the average fitting method for the calculated acceleration response spectrum.The design ground motion acceleration response spectrum of the project site is taken as follows:(9)SaT=AmaxβT,αmax=Amaxβmaxg,where Amax is the design peak ground acceleration, αmax the maximum value of seismic influence coefficient, g is the gravitational acceleration, and βT is given in the form of《Code for seismic design of Highway Engineering》(JTG B02-2013):(10)βT=1+βmax−1TT1,0≤T<T1s,βmax,T1≤T≤Tgs,βmaxTg/Tγ,T>Tgs,where T is the natural vibration period of the structure, βmax is the maximum value of the response spectrum, T1 is the initial period of the platform section with the maximum response spectrum, Tg is the characteristic period, and γ is the attenuation index of the falling section of the response spectrum.


### 4.4. Design Acceleration Time History

In this work, the artificial ground motion is synthesized according to the calculated peak value and response spectrum results. Formula (10) is used to calculate the horizontal ground motion acceleration response spectrum results combined with the 50-year exceedance probability of the project site of 10%. The relative error between the time history response spectrum and the target spectrum is less than 5%. The synthesis takes 0.02 seconds as the step length, and the envelope function mainly considers the influence of long-period ground motion. By calculating the equivalent magnitude and equivalent distance of ground motion corresponding to the natural vibration period of the structure, the envelope parameters used are determined. The horizontal acceleration time history curve of each exceedance probability and its fitting with the target spectrum are shown in Figures 6 and 7.Figure 6
Time history of horizontal design ground motion acceleration of 10% section II foundation surface.Figure 7
Fitting of the response spectrum.


## 4.1. Synthetic Ground Motion Time History of Bedrock

In the seismic response time history analysis, the selection of an appropriate ground motion acceleration time process is very important. Although the number of actual seismic records has greatly increased in the past few decades, the site conditions of the recording sites may be very different from the construction site conditions we are concerned about. We need a set of ground motion samples that meet the same statistical characteristics as the seismic response input. In this section, based on the above-ground motion characteristic parameters of the project site, the triangular series superposition simulation method is used to synthesize the bedrock ground motion time history. The calculation model is as follows:(1)x″t=ft∑k=0nCkcosωkt+φk=ftat,(2)at=∑k=0nCkcosωkt+φk.In the formula,φk is the random phase angle uniformly distributed in the 0,2π interval; Ck and ωk is the amplitude and frequency of the k-th frequency component, respectively; ft is the strength envelope function, which is a definite function of time; and at is a stationary Gaussian process.It can be seen from equation (1) that the synthetic acceleration time history x″t with the characteristics of a nonstationary random process is the product of a Gaussian process at with the characteristics of a stationary random process and an intensity envelope function ft.The coefficientCk in equation (2) can be determined by a given power spectral density function Sωk , where we obtain as follows:(3)Ck=4SωkΔω,ωk=2πkT,Δω=2πT,where t is the total holding time of the stationary random process. In order to use the acceleration response spectrum as the target spectrum of artificial acceleration time history, the following approximate relationship between the response spectrum and power spectrum can be used to replace Sωk in equation (3);(4)Sωk=−ζπωSaTω21ln−π/ωTln1−P.In the formula,SaTω is the given target acceleration response spectrum, ζ is the damping ratio, t is the duration, and P is the exceedance probability of the response.The intensity envelope functionft is related to the magnitude, epicenter distance, geological structure background, site soil conditions, and other factors and reflects the unstable characteristics of ground motion with time. Generally, the strength envelope function ft is expressed by the following piecewise function:(5)ft=t/t12,0≤t≤t1,1,t1<t≤t2,e−ct−t2,t2<t≤td,where 0,t1 is the rising section, t1,t2 is the stationary section, t2,td is the attenuation section, and C is the attenuation coefficient. Since the relationship between the response spectrum and the power spectrum represented by equation (4) is approximate, the response spectrum calculated according to the initial time history is generally only approximate to the target response spectrum. In order to improve the fitting accuracy, iterative adjustment is required. When the difference between the calculated response spectrum and the target response spectrum is less than the control accuracy (take 5%), the iteration is stopped.


## 4.2. Dynamic Response Analysis of Site Soil

### 4.2.1. Calculation Model of Dynamic Response of the Site Soil Layer

The results of site engineering seismic conditions show that within the scope of the project site, the changes in medium characteristics and terrain are not very significant along the horizontal direction. Therefore, the influence of site conditions on seismic ground motion is considered based on a one-dimensional site model. The equivalent linearization method of one-dimensional soil shear dynamic response analysis is used in the study of the influence of the one-dimensional site model on earthquake ground motion. Its basic principle is as follows:It is assumed that the shear wave is incident vertically from the viscoelastic semi-infinite bedrock space into the horizontal layered (N-layer) nonlinear soil and propagates upward. For this calculation model, according to the wave propagation theory, the dynamic response value of the field ground medium can be calculated by using the time-frequency transformation technology, combined with the complex damping simulation of the nonlinear characteristics of soil and the equivalent linearization method.A shear harmonic is set to incident vertically upward from the calculated base and propagate in the soil layer. According to the wave theory and complex damping theory, the medium motion in each soil layer must meet the wave equation:(6)ρj∂2Ujx.t∂t2=Gjc∂2Ujx,t∂x2,where Ujx,t is the displacement value of medium reaction in the soil layer j, ρj is the density of the medium in the soil layer j, Gjc is the dynamic complex shear modulus of the medium in the soil layer j, Gjc is given by the following equation:(7)Gjc=1+2λjγjeiGjdγjeGjo,whereGjo is the maximum dynamic shear modulus of medium in soil layer, jGjdγje,λjγje is the dimensionless coefficient of equivalent dynamic shear modulus and hysteretic damping ratio of the medium in soil layer, and jγje is the maximum dynamic shear modulus of the medium in soil layer jThe medium motion between soil layers meets the conditions of displacement continuity and stress continuity.(8)Ujx,tx=Hj=Uj+1x,tx=0,τjx,tx=Hj=τj+1x,tx=0,τ1x,tx=0=0.In the formula,Hj is the layer thickness of the j-th soil layer, and it is specified that the vertically downward direction is the positive direction of the X coordinate and the coordinate origin is placed at the top of each soil layer. Through formula (6), the frequency domain value of the medium reaction quantity in the soil layer can be obtained by using the known calculated base incident wave value, and then, the time domain value of the medium reaction quantity in the soil layer can be obtained by using the Fourier transform method.Considering the nonlinear characteristics of soil, the dimensionless coefficient of equivalent dynamic shear modulus and hysteretic damping ratio of each soil layer are functions of equivalent shear strain. Therefore, in the actual calculation, the initial equivalent dynamic shear strain of the medium reaction in each soil layer is assumed and the response calculation is carried out by using the abovementioned method. Then, the maximum shear strain response of the medium at the midpoint of each soil layer is calculated. Finally, the maximum shear strain of the medium reaction at the midpoint of each soil layer multiplied by the reduction coefficient (0.65) is taken as the calculated value of the equivalent shear strain of the medium in the soil layer. We compare the equivalent shear strain used in the calculation with the equivalent dynamic shear modulus and hysteretic damping ratio corresponding to the calculated equivalent shear strain. If the relative error is less than the given allowable error (0.05), it is considered that the consideration of the nonlinear characteristics of soil meets the requirements. Otherwise, we replace the initial equivalent shear strain value with the latest calculated equivalent shear strain value and repeat the abovementioned calculation process until the relative error is less than the allowable error.


### 4.2.2. Determination of Dynamic Parameters of the Site Soil Layer

Focal characteristics, propagation path, and local site conditions are the main factors affecting the ground motion of the site, and the influence of local site conditions is particularly prominent in a small area. Therefore, the engineering geological exploration in the proposed site and the determination of the type, stratification, and thickness of the site soil are essential work in the seismic response calculation. The analysis and determination of dynamic mechanical indexes of site soil (such as bulk density, shear wave velocity, shear modulus ratio curve, and damping ratio curve) are very important for site dynamic response analysis.In this work, the existing data are sorted out according to the requirements of soil seismic response analysis and calculation. According to the borehole data of the site, a borehole with shear wave velocity data is taken as the calculation control point for each station site and station-to-station interval. In the seismic response calculation, the medium with a shear wave velocity greater than 500 m/s is determined as the seismic input interface and the shear wave velocity of the soil below this layer is greater than 500 m/s.


## 4.2.1. Calculation Model of Dynamic Response of the Site Soil Layer

The results of site engineering seismic conditions show that within the scope of the project site, the changes in medium characteristics and terrain are not very significant along the horizontal direction. Therefore, the influence of site conditions on seismic ground motion is considered based on a one-dimensional site model. The equivalent linearization method of one-dimensional soil shear dynamic response analysis is used in the study of the influence of the one-dimensional site model on earthquake ground motion. Its basic principle is as follows:It is assumed that the shear wave is incident vertically from the viscoelastic semi-infinite bedrock space into the horizontal layered (N-layer) nonlinear soil and propagates upward. For this calculation model, according to the wave propagation theory, the dynamic response value of the field ground medium can be calculated by using the time-frequency transformation technology, combined with the complex damping simulation of the nonlinear characteristics of soil and the equivalent linearization method.A shear harmonic is set to incident vertically upward from the calculated base and propagate in the soil layer. According to the wave theory and complex damping theory, the medium motion in each soil layer must meet the wave equation:(6)ρj∂2Ujx.t∂t2=Gjc∂2Ujx,t∂x2,where Ujx,t is the displacement value of medium reaction in the soil layer j, ρj is the density of the medium in the soil layer j, Gjc is the dynamic complex shear modulus of the medium in the soil layer j, Gjc is given by the following equation:(7)Gjc=1+2λjγjeiGjdγjeGjo,whereGjo is the maximum dynamic shear modulus of medium in soil layer, jGjdγje,λjγje is the dimensionless coefficient of equivalent dynamic shear modulus and hysteretic damping ratio of the medium in soil layer, and jγje is the maximum dynamic shear modulus of the medium in soil layer jThe medium motion between soil layers meets the conditions of displacement continuity and stress continuity.(8)Ujx,tx=Hj=Uj+1x,tx=0,τjx,tx=Hj=τj+1x,tx=0,τ1x,tx=0=0.In the formula,Hj is the layer thickness of the j-th soil layer, and it is specified that the vertically downward direction is the positive direction of the X coordinate and the coordinate origin is placed at the top of each soil layer. Through formula (6), the frequency domain value of the medium reaction quantity in the soil layer can be obtained by using the known calculated base incident wave value, and then, the time domain value of the medium reaction quantity in the soil layer can be obtained by using the Fourier transform method.Considering the nonlinear characteristics of soil, the dimensionless coefficient of equivalent dynamic shear modulus and hysteretic damping ratio of each soil layer are functions of equivalent shear strain. Therefore, in the actual calculation, the initial equivalent dynamic shear strain of the medium reaction in each soil layer is assumed and the response calculation is carried out by using the abovementioned method. Then, the maximum shear strain response of the medium at the midpoint of each soil layer is calculated. Finally, the maximum shear strain of the medium reaction at the midpoint of each soil layer multiplied by the reduction coefficient (0.65) is taken as the calculated value of the equivalent shear strain of the medium in the soil layer. We compare the equivalent shear strain used in the calculation with the equivalent dynamic shear modulus and hysteretic damping ratio corresponding to the calculated equivalent shear strain. If the relative error is less than the given allowable error (0.05), it is considered that the consideration of the nonlinear characteristics of soil meets the requirements. Otherwise, we replace the initial equivalent shear strain value with the latest calculated equivalent shear strain value and repeat the abovementioned calculation process until the relative error is less than the allowable error.


## 4.2.2. Determination of Dynamic Parameters of the Site Soil Layer

Focal characteristics, propagation path, and local site conditions are the main factors affecting the ground motion of the site, and the influence of local site conditions is particularly prominent in a small area. Therefore, the engineering geological exploration in the proposed site and the determination of the type, stratification, and thickness of the site soil are essential work in the seismic response calculation. The analysis and determination of dynamic mechanical indexes of site soil (such as bulk density, shear wave velocity, shear modulus ratio curve, and damping ratio curve) are very important for site dynamic response analysis.In this work, the existing data are sorted out according to the requirements of soil seismic response analysis and calculation. According to the borehole data of the site, a borehole with shear wave velocity data is taken as the calculation control point for each station site and station-to-station interval. In the seismic response calculation, the medium with a shear wave velocity greater than 500 m/s is determined as the seismic input interface and the shear wave velocity of the soil below this layer is greater than 500 m/s.


## 4.3. Site Design Ground Motion Parameters

On the basis of the seismic response calculation results of the site soil layer obtained in the previous section, the design ground motion parameters of the project site in this section will be given, including the design peak ground motion acceleration and acceleration response spectrum. The design peak ground acceleration is determined by comprehensively considering the peak value of the sample and the short-period acceleration response spectrum. The design acceleration response spectrum is determined by using the average fitting method for the calculated acceleration response spectrum.The design ground motion acceleration response spectrum of the project site is taken as follows:(9)SaT=AmaxβT,αmax=Amaxβmaxg,where Amax is the design peak ground acceleration, αmax the maximum value of seismic influence coefficient, g is the gravitational acceleration, and βT is given in the form of《Code for seismic design of Highway Engineering》(JTG B02-2013):(10)βT=1+βmax−1TT1,0≤T<T1s,βmax,T1≤T≤Tgs,βmaxTg/Tγ,T>Tgs,where T is the natural vibration period of the structure, βmax is the maximum value of the response spectrum, T1 is the initial period of the platform section with the maximum response spectrum, Tg is the characteristic period, and γ is the attenuation index of the falling section of the response spectrum.


## 4.4. Design Acceleration Time History

In this work, the artificial ground motion is synthesized according to the calculated peak value and response spectrum results. Formula (10) is used to calculate the horizontal ground motion acceleration response spectrum results combined with the 50-year exceedance probability of the project site of 10%. The relative error between the time history response spectrum and the target spectrum is less than 5%. The synthesis takes 0.02 seconds as the step length, and the envelope function mainly considers the influence of long-period ground motion. By calculating the equivalent magnitude and equivalent distance of ground motion corresponding to the natural vibration period of the structure, the envelope parameters used are determined. The horizontal acceleration time history curve of each exceedance probability and its fitting with the target spectrum are shown in Figures 6 and 7.Figure 6
Time history of horizontal design ground motion acceleration of 10% section II foundation surface.Figure 7
Fitting of the response spectrum.


## 5. Response Analysis of the Structure

### 5.1. Finite Element Model

The finite element model of the subway station is established by Midas-GTS software, and the response analysis of the subway station under seismic load is carried out using the time history method. The dimension of the numerical model is selected as 100 m × 40 m × 47 m to include most of the soil affected by the excavation. That is, the side artificial boundary of the calculation model is 3 times the horizontal effective width of the subway station. The lower boundary is to the equivalent bedrock surface, and the upper surface is taken to the actual surface. Horizontal and vertical seismic acceleration loads are applied at the bottom of the model, and the peak value of vertical acceleration is 65% of the horizontal acceleration. Free field boundary conditions are set around the calculation model. To prevent seismic wave reflection, the unit width of the free boundary is set to 1000 km and the bottom of the model is a fixed constraint. To ensure calculation accuracy and reduce the calculation time as much as possible, the grid near the main structure of the station is densified. The stratum is simulated by a three-dimensional solid element, the station secondary lining is simulated by a plate element, and the middle column and longitudinal beam are simulated by a beam element. The Mohr–Coulomb model is adopted for the soil constitutive relationship and the elastic-plastic model is adopted for station lining. The calculation model is shown in Figure8.Figure 8
Calculation model.


### 5.2. Result Analysis

The safety of the station under earthquake is directly affected by the relative displacement of the top and bottom plates of the subway station. It can be seen from Figure9 that under the seismic conditions of fortification intensity, the maximum horizontal displacement of the bottom plate is 18.8 mm. The horizontal displacement of the middle plate and the bottom plate are 16.7 and 15.0 mm, respectively. Therefore, the interlayer displacement between the top plate and the middle plate is 2.1 mm and the interlayer displacement angle is 1/2650. The interlayer displacement between the middle baseplates is 1.7 mm and the interlayer displacement angle is 1/4117, see Figure 10. Under the action of a fortification earthquake, the inter story displacement angle of Taiyuan Street station is less than 1/550, indicating that the structure meets the seismic requirements.Figure 9
Horizontal displacement at the moment of maximum relative deformation of the top and bottom plates. (a) Overall horizontal displacement. (b) Horizontal displacement of the structure.
(a)(b)Figure 10
Relative displacement history curve of the top and bottom plates.


## 5.1. Finite Element Model

The finite element model of the subway station is established by Midas-GTS software, and the response analysis of the subway station under seismic load is carried out using the time history method. The dimension of the numerical model is selected as 100 m × 40 m × 47 m to include most of the soil affected by the excavation. That is, the side artificial boundary of the calculation model is 3 times the horizontal effective width of the subway station. The lower boundary is to the equivalent bedrock surface, and the upper surface is taken to the actual surface. Horizontal and vertical seismic acceleration loads are applied at the bottom of the model, and the peak value of vertical acceleration is 65% of the horizontal acceleration. Free field boundary conditions are set around the calculation model. To prevent seismic wave reflection, the unit width of the free boundary is set to 1000 km and the bottom of the model is a fixed constraint. To ensure calculation accuracy and reduce the calculation time as much as possible, the grid near the main structure of the station is densified. The stratum is simulated by a three-dimensional solid element, the station secondary lining is simulated by a plate element, and the middle column and longitudinal beam are simulated by a beam element. The Mohr–Coulomb model is adopted for the soil constitutive relationship and the elastic-plastic model is adopted for station lining. The calculation model is shown in Figure8.Figure 8
Calculation model.


## 5.2. Result Analysis

The safety of the station under earthquake is directly affected by the relative displacement of the top and bottom plates of the subway station. It can be seen from Figure9 that under the seismic conditions of fortification intensity, the maximum horizontal displacement of the bottom plate is 18.8 mm. The horizontal displacement of the middle plate and the bottom plate are 16.7 and 15.0 mm, respectively. Therefore, the interlayer displacement between the top plate and the middle plate is 2.1 mm and the interlayer displacement angle is 1/2650. The interlayer displacement between the middle baseplates is 1.7 mm and the interlayer displacement angle is 1/4117, see Figure 10. Under the action of a fortification earthquake, the inter story displacement angle of Taiyuan Street station is less than 1/550, indicating that the structure meets the seismic requirements.Figure 9
Horizontal displacement at the moment of maximum relative deformation of the top and bottom plates. (a) Overall horizontal displacement. (b) Horizontal displacement of the structure.
(a)(b)Figure 10
Relative displacement history curve of the top and bottom plates.


## 6. Conclusions

(1)
The dynamic elastic modulus of sand decreases with the increase of dynamic elastic strain and increases with the increase of confining pressure. The influence of confining pressure on the damping ratio of sand is slight, but it can still be seen that the damping ratio decreases with the increase of confining pressure, especially when the dynamic elastic strain is relatively low.(2)
The interlayer displacement between the top and middle plate is 2.1 mm and the interlayer displacement angle is 1/2650. The interlayer displacement between the middle and bottom plate is 1.7 mm and the interlayer displacement angle is 1/4117. Under the action of a fortification earthquake, the interlayer displacement angles of Taiyuan Street station are less than 1/550, indicating that the structure meets the seismic requirements.



---

*Source: 5294597-2022-09-24.xml*