Numerical Simulation and In-Situ Measurement of Ground-Borne Vibration Due to Subway System

Jinping Yang 1, Peizhen Li 1,2 ID and Zheng Lu 1,2,* ID 1 Research Institute of Structural Engineering and Disaster Reduction, Tongji University, Shanghai 200092, China; 1410231@tongji.edu.cn (J.Y.); lipeizh@tongji.edu.cn (P.L.) 2 State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China * Correspondence: luzheng111@tongji.edu.cn; Tel.: +86-216-598-6186


Introduction
In recent decades, with the rapid development of urban modernization, the impact of traffic system induced vibration on urban life and working environment has been brought to public attention since vibration has been listed as one of the seven major environmental hazards in the world [1].Subways, key parts of the urban traffic, are highly demanding on the environment, particularly in major cities, e.g., Shanghai and Beijing, China.Therefore, the vibration induced by the passage of underground trains is a major environmental concern in urban areas [2].
The influences of subway induced vibration is mainly reflected in the following areas: the influence of vibration on the surrounding buildings, especially historical architecture, such as the collapse of the old church caused by vibration.Vibration may also interfere with sensitive equipment used in scientific research and high-tech industries [3].Moreover, vibration and secondary noise affect people's work, life and even their health.Whole-body human exposure to vibration indoors has been evaluated in the frequency range 1-80 Hz [4], which is consistent with the range of subway induced vibration [2].Therefore, it is imperative to focus on the subway induced vibration.
Great efforts have been made on the measurements of subway induced vibration to study the vibration characteristics, the vibration causes, the propagation path and the controlling means [5].Vibration measurements are performed on the building foundation slab and in open fields adjacent to the building [3], the subway train-induced vibrations [6], viaduct track and nearby ground [7].Vibration measurement systems are various, including a remote system for the online monitoring [8], real-time monitoring of adjacent cross operation subway [9], and field monitoring at permafrost testing section [10].Moreover, effective measures to mitigate vibration are developed and suggested based on field or test results, such as a floating stab tracks [11], rail pads [12], and wave barrier backfill material [13].Besides, track vibration isolation [14] and energy-absorbing rubber [15,16] are effective strategies to mitigate the rail-way infrastructure issues, and structural resonant behavior can be avoided [17][18][19].Moreover, structural vibration control methods [20][21][22], such as tuned mass damper [23][24][25], nonlinear energy sink [26,27], particle damper [28][29][30], etc., may also be applied to suppress subway induced vibrations.Some interrelated research, theoretical studies and numerical analysis of subway induced vibration have recently received attention [31].Numerical modeling includes: two dimensional (2D) models [32], two-and-a-half-dimensional (2.5D) models [33] and three-dimensional (3D) models [34].2D analytical methods have the advantage of high computational efficiency and are suitable for simple geometries.2.5D models have received wide concern and many studies have been performed, including on the 2.5D finite/infinite element approach [35] and the 2.5D finite element coupled boundary element BEM model [36].Other novel numerical methods, such as a closed-form semi-analytical solution in a half-space [37], have also been developed for the prediction of subway induced vibrations.3D models have high accuracy, however, are more computationally expensive [38].
It has proven effective to combine numerical studies together with experimental results in studies on the subway induced vibration.Chiacchiari [39] discussed measurement methods for rail corrugation and proposed an application method for rail irregularities.Kaynia [40] performed a measurement of the vibration induced railway embankment by high speed trains and developed a rigorous numerical model for the prediction of ground response.
It can be seen from the above studies that the issue of subway induced vibration is a sustainable topic, causing annoyance and adverse physical, physiological, and psychological effects on humans in dense urban environments [5].Consequently, the adopted approaches and methods in the vibration mitigation must be sustainable to reduce damage and adverse effects on humans, and it is of great significance to focuses on methods to enhance the sustainability in design of newly built metro depots.
Based on the basic principles of the infinite element method, a coupled finite element-infinite element boundary method that combines the finite element method in the key areas and infinite element method in the far field is introduced in Section 2 to address the issue of subway induced vibration.Then, the finite element method performed by the harmonic analysis based on the commercial software ABAQUS is verified by comparing with the thin layer method.An interface program is developed to realize the frequency domain calculation on the MATLAB platform and it is verified by the transfer function method.In Section 3, a field measurement of practical engineering is performed on a track bed-tunnel lining-surrounding formation located on Line 2, Shanghai Metro, and the corresponding simulation model is established based on the coupled method.The numerical results of vibration level and spectrum characteristic, obtained by the developed interface programs, are discussed.Then, a simplified prediction formula of ground vibration level is proposed based on the measurement results.Finally, parametric study is conducted in an elastic half-space simulation to investigate the influence of model widths and depths with the combined method on the numerical results in Section 4.

Basic Principles and Realization of Infinite Element Method in Programs
Infinite element, a special kind of finite element, is a unit that tends to infinity geometrically.It is an effective supplement to the finite element in dealing with the unbounded domain issues, and can be seamlessly connected to the finite element.
The directionality is the big difference between infinite and finite units and the two sides cannot intersect in the infinite direction.Therefore, the infinite element must be divergent in the infinite direction mathematically.Only under these conditions can the quantities be monotonic in the process of tending towards infinity.In the dynamic analysis of infinite large medium (e.g., foundations), artificial boundary is affected by the wave reflection, and therefore the energy is transmitted back to the analytical zone and cannot be propagated to infinity.Damping force is introduced to absorb the radiant energy of planar body wave in the infinite boundary.
The damping stress is introduced to the boundary of x = L and can be expressed as follows: where .u x , .u y and .
u z are the vibration velocities in x, y and z direction, respectively.The constants d p and d s , proportional to the velocity, are selected to avoid the energy reflection of longitudinal and shear waves to the infinite zones.They are generally expressed as: where c p and c s are the propagation velocities of longitudinal wave and transverse wave, respectively.The infinite elements include plain strain element, plain stress element, axial symmetry element and three-dimensional infinite element in ABAQUS program.The first surface of the infinite element is required to be the interface between finite element and infinite element to ensure the direction of the infinite element extending from near field to far field.In addition, the infinite elements cannot be intersected with each other in the extension direction.
The coupled finite-infinite-infinite element boundary is performed in the numerical analysis.The focused areas of interest (near field) are simulated by finite elements and the far field by infinite elements [41].Finite elements with certain nodes are firstly defined, and the unit type is then changed in INP file to accomplish the definition of infinite elements, which cannot be defined directly in ABAQUS program.The simulation model is conducted in the harmonic analysis based on the commercial software ABAQUS.

Verification of Infinite Element Method in Programs
It is an effective way to test the applicability of the presented model in different methods.The proposed finite element coupled infinite element boundary method in ABAQUS software is compared with the thin layer method (TLM), a semi-analytical and semi-numerical approach for the analysis of elastic wave propagation in stratified media [42].
An elastic half-space simulation is taken as an example in the subway induced vibration models.The basic parameters of the numerical model are as follows: the mass density ρ is 2000 kg/m 3 ; shear wave velocity Cs is 500 m/s; Poisson's ratio µ is 0.4; and the material damping ratio ξ of soil is 0.05.The depth of soil layer is 80 m and the base is assumed as fixed rock.The two-dimensional numerical model is built and half-structure is utilized based on symmetry in the present study.The cosine function F = F 0 cos(2π f t), with the amplitude F 0 5 kN, is applied to the models, and the applied frequency range f is from 1 Hz to 80 Hz according to the required frequency band under the subway-induced vibration.The ground evaluation points are within the range of 2-30 m from the distance of load application with total 15 points and 2 m interval and the vertical vibrating displacement of these points are obtained and discussed in this section.
The base of the numerical model is fixed and the horizontal degrees are constrained on the right side.The model can be expressed as L_H, where L denotes the width of the model and H is the height (depth).The model 1280_80 notes the model width is 1280 m and the depth is 80 m.ABAQUS is abbreviated to ABA in the numerical figures.
Figure 1 presents the comparison attenuation curves at the typical frequencies of TLM and ABAQUS results with element dimension 1 m, which can be drawn that the numerical results of the two methods are similar under the low frequency load, such as 8 Hz.There are differences under the high frequency load due to the larger element dimensions.However, the differences become smaller with the finer element dimensions [35,41] in ABAQUS software, as shown in Figure 2, indicating that the finite element method is verified applicable in ABAQUS software.
Sustainability 2018, 10, x FOR PEER REVIEW 4 of 20 frequency range f is from 1 Hz to 80 Hz according to the required frequency band under the subway-induced vibration.The ground evaluation points are within the range of 2-30 m from the distance of load application with total 15 points and 2 m interval and the vertical vibrating displacement of these points are obtained and discussed in this section.The base of the numerical model is fixed and the horizontal degrees are constrained on the right side.The model can be expressed as L_H, where L denotes the width of the model and H is the height (depth).The model 1280_80 notes the model width is 1280 m and the depth is 80 m.ABAQUS is abbreviated to ABA in the numerical figures.
Figure 1 presents the comparison attenuation curves at the typical frequencies of TLM and ABAQUS results with element dimension 1 m, which can be drawn that the numerical results of the two methods are similar under the low frequency load, such as 8 Hz.There are differences under the high frequency load due to the larger element dimensions.However, the differences become smaller with the finer element dimensions [35,41] in ABAQUS software, as shown in Figure 2, indicating that the finite element method is verified applicable in ABAQUS software.

The Realization of Infinite Element Method in the Frequency Domain
To make the simulation results statistically significant, many excitations are applied to the twodimensional, or more complex, models in the subway induced vibration, since the subway incentive is a random process.The required calculation process becomes time-consuming in the time domain analysis.Moreover, the frequency characteristics of the vibration is important in the subway induced vibration.Therefore, the numerical analysis conducted in the frequency domain is of great frequency range f is from 1 Hz to 80 Hz according to the required frequency band under the subway-induced vibration.The ground evaluation points are within the range of 2-30 m from the distance of load application with total 15 points and 2 m interval and the vertical vibrating displacement of these points are obtained and discussed in this section.The base of the numerical model is fixed and the horizontal degrees are constrained on the right side.The model can be expressed as L_H, where L denotes the width of the model and H is the height (depth).The model 1280_80 notes the model width is 1280 m and the depth is 80 m.ABAQUS is abbreviated to ABA in the numerical figures.
Figure 1 presents the comparison attenuation curves at the typical frequencies of TLM and ABAQUS results with element dimension 1 m, which can be drawn that the numerical results of the two methods are similar under the low frequency load, such as 8 Hz.There are differences under the high frequency load due to the larger element dimensions.However, the differences become smaller with the finer element dimensions [35,41] in ABAQUS software, as shown in Figure 2, indicating that the finite element method is verified applicable in ABAQUS software.

The Realization of Infinite Element Method in the Frequency Domain
To make the simulation results statistically significant, many excitations are applied to the twodimensional, or more complex, models in the subway induced vibration, since the subway incentive is a random process.The required calculation process becomes time-consuming in the time domain analysis.Moreover, the frequency characteristics of the vibration is important in the subway induced vibration.Therefore, the numerical analysis conducted in the frequency domain is of great

The Realization of Infinite Element Method in the Frequency Domain
To make the simulation results statistically significant, many excitations are applied to the two-dimensional, or more complex, models in the subway induced vibration, since the subway incentive is a random process.The required calculation process becomes time-consuming in the time domain analysis.Moreover, the frequency characteristics of the vibration is important in the subway induced vibration.Therefore, the numerical analysis conducted in the frequency domain is of great significance.The frequency analysis is a process of data operation, rather than the step-by-step solution process in the time history analysis.The calculation time is then saved effectively.
In this study, an interface program is developed in the commercial software ABAQUS environment.It can automatically read the simulation result files (dat file) in the ABAQUS harmonic analysis, and then the frequency domain calculation is realized by MATLAB.The main steps of frequency domain analysis could be summarized by the following: Step 1.The coupled finite element-infinite element boundary simulation model is built based on the commercial software ABAQUS.
Step 2. The harmonic analyses are conducted based on the commercial software ABAQUS to obtain the amplitude and phase data of the evaluation measures on the ground in the subway induced vibration.Step 3. The numerical results are saved as dat file type in ABAQUS, and then read automatically by the developed interface programs to MATLAB.Step 4. Frequency domain analysis is performed based on MATLAB platform and the responses, such as acceleration and displacement time-histories of the evaluation points are obtained.
The simulation results acquired in Steps 1 and 2 are verified by TLM in Section 2.2.The time history with certain frequency is obtained based on the amplitude and phase, which are calculated in Step 2.Then, the total responses of evaluation points with each frequency are obtained based on the superposition principle.The frequency domain analysis in Step 4 is verified by the transfer function method, whose principle is expressed as follows: If the transfer function H(ω) is known, then the responses of the system can be calculated by Equations ( 4) and (5).
where X(ω) is the frequency response functions and Y(ω) is the output Fourier transform.Therefore, the responses of system can be calculated under any excitation X(t), provided the transfer function H(ω) is known.
The correctness and the application of the program is verified by a numerical simulation, a track bed-tunnel lining-surrounding formation model, as shown in Figure 3a.Displacement load, as expressed in Equation (6), is applied to the center of the tunnel track-bed (A point) in the vertical direction.The applied displacement curve is shown in Figure 3b, and the displacement time-history of the ground surface B point is required to be calculated.Table 1 lists the soil parameters of different soil layers in this model.
Figure 4 shows the comparison displacement time history of Point B by the developed interface program in frequency domain analysis and the transfer function method.The two curves almost coincide with each other, indicating that the frequency analysis is applicable in the interface programs.There is little difference between the two results in the amplified graph, as shown in Figure 4b.The reason is that the sampling points are limited in the frequency analysis, resulting in frequency leakage and numerical calculation errors.

Vibration Measurement
Field measurements were performed on a track bed-tunnel lining-surrounding formation in a certain park project of Line 2, Shanghai Metro, which is an ideal test site due to the few surroundings and buildings to study the vibration attenuation characteristics of the ground surface under the subway induced vibration.

Vibration Measurement
Field measurements were performed on a track bed-tunnel lining-surrounding formation in a certain park project of Line 2, Shanghai Metro, which is an ideal test site due to the few surroundings and buildings to study the vibration attenuation characteristics of the ground surface under the subway induced vibration.The acceleration time-histories of the central tunnel and the typical ground points were measured and acquired in this test.The track bed-tunnel lining-surrounding formation model and the arrangement of measuring points in lateral direction, which were 0 m, 5 m, 11 m, 15 m, 25 m and 40 m from the central tunnel, are shown in Figure 5.
The outer radius of the tunnel is 3.1 m and inner radius is 2.75 m.The tunnel depth, the distance from the top of the tunnel to the surface of the tunnel, is 13.8 m.The elastic modulus, the density and Poisson's ratio of concrete are 34,500 MPa, 2500 kg/m 3 and 0.25, respectively.The typical clay in Shanghai, not homogeneous in the depth direction, was employed in this field test and numerical simulation.The soil parameters of different soil layers are shown in Table 1  Poisson's ratio of concrete are 34,500 MPa, 2500 kg/m 3 and 0.25, respectively.The typical clay in Shanghai, not homogeneous in the depth direction, was employed in this field test and numerical simulation.The soil parameters of different soil layers are shown in Table 1 in Section 2.3.The 31 accelerations of the evaluation points in the center of the track bed were measured and acquired in this field test, with total excitation time 10 s and sample frequency 200 Hz.The 31 accelerations correspond to 31 cases.The typical acceleration time history of the track bed, its corresponding Fourier spectra and 1/3 octave frequency are provided in Figure 6a-c, in which it can be observed that the dominant frequencies of the tunnel track bed induced by the subway are above 40 Hz.The measured results agree with other tests [43] that the vibration frequency induced by the subway is high.Vibration level of each recorded accelerations is shown in Figure 6d based on the Shanghai regulation, in which it can be seen that the vibration level is within 80-85 dB and the mean value is 82.4 dB.

The Coupled Finite Element-Infinite Element Boundary Model
In this engineering, the minimum wave velocity of the soil is 100 m/s, the width of the finite element region is designed as 100 m, and the model depth is 100 m.Infinite element boundary is applied to the right side and the bottom of the models.The element size in this engineering along the horizontal direction and vertical direction within 50 m depth is designed as 0.2 m, while other regions are meshed with 0.4 m.The tunnel is analyzed with plane strain solid elements, with the element size 0.1 m × 0.1 m (width and height).The measured 31 accelerations of the evaluation points in the center of the track bed are applied vertically to this position in the simulation model.The simulation frequency is 0-100 Hz, with 0.1 Hz interval.The simulation points are within 0-40 m of the center of the tunnel, coincident to the measure points displayed in Figure 5.The meshed model of track bed-tunnel lining-surrounding formation is shown in Figure 7.The numerical results are obtained by the harmonic analysis in ABAQUS software and frequency analysis in MATLAB, connecting by the developed interface programs, as introduced in Section 2.3.Soil behaviors are simplified in this simulation model and soil properties are based on the practical engineering, as shown in Table 1.

The Coupled Finite Element-Infinite Element Boundary Model
In this engineering, the minimum wave velocity of the soil is 100 m/s, the width of the finite element region is designed as 100 m, and the model depth is 100 m.Infinite element boundary is applied to the right side and the bottom of the models.The element size in this engineering along the horizontal direction and vertical direction within 50 m depth is designed as 0.2 m, while other regions are meshed with 0.4 m.The tunnel is analyzed with plane strain solid elements, with the element size 0.1 m × 0.1 m (width and height).The measured 31 accelerations of the evaluation points in the center of the track bed are applied vertically to this position in the simulation model.The simulation frequency is 0-100 Hz, with 0.1 Hz interval.The simulation points are within 0-40 m of the center of the tunnel, coincident to the measure points displayed in Figure 5.The meshed model of track bedtunnel lining-surrounding formation is shown in Figure 7.The numerical results are obtained by the harmonic analysis in ABAQUS software and frequency analysis in MATLAB, connecting by the developed interface programs, as introduced in Section 2.3.Soil behaviors are simplified in this simulation model and soil properties are based on the practical engineering, as shown in Table 1.

Numerical Results
The simulation attenuation of peak vibration levels in the vertical direction with respect to distance under the 31 measured accelerations excitations are shown in Figure 8a,b.In addition, the mean attenuations of peak vibration levels, which are obtained by the measurement and simulation, are compared in Figure 8c.The following conclusions can be drawn.
The attenuation curves of vibration levels generally decrease with the increasing distance away from the central tunnel.However, there exists an amplifying area of vibration that the vibration level

Numerical Results
The simulation attenuation of peak vibration levels in the vertical direction with respect to distance under the 31 measured accelerations excitations are shown in Figure 8a,b.In addition, the mean attenuations of peak vibration levels, which are obtained by the measurement and simulation, are compared in Figure 8c.The following conclusions can be drawn.
The attenuation curves of vibration levels generally decrease with the increasing distance away from the central tunnel.However, there exists an amplifying area of vibration that the vibration level at 10 m from the center is about 5 dB greater than that at 15 m.This amplification effect due to the wave propagation and reflection is consistent with the results of measurement and numerical simulation obtained by other researchers [44].If the reflection frequency is dominant and obvious in the applied frequency components, the overall vibration will reflect due to the wave superposition in this area.A part of the vibration wave is directly transmitted through the bound site, and part of the wave first spreads to the ground surface and then propagates from the ground to the lateral sides, forming a superposition effect in this area.This may result in a larger vibration in this field than that directly above the tunnel.The wave propagation phenomenon is clearly observed in the measured ground data by Takemiya [7].
Moreover, the amplification phenomenon is obvious in some cases, such as in Case 11, while some are not apparent, such as in Case 29, as shown in Figure 8b, indicating that the amplification has relationship with the spectrum of the applied excitation accelerations.
The numerical mean value of vibration level basically maintains within the measured range, which is a little smaller than the mean value of experimental results, as shown in Figure 8c.The amplification of vibration level is also obtained in the simulation, about 15 m away from the center of the tunnel.
Sustainability 2018, 10, x FOR PEER REVIEW 9 of 20 larger vibration in this field than that directly above the tunnel.The wave propagation phenomenon is clearly observed in the measured ground data by Takemiya [7].Moreover, the amplification phenomenon is obvious in some cases, such as in Case 11, while some are not apparent, such as in Case 29, as shown in Figure 8b, indicating that the amplification has relationship with the spectrum of the applied excitation accelerations.
The numerical mean value of vibration level basically maintains within the measured range, which is a little smaller than the mean value of experimental results, as shown in Figure 8c.The amplification of vibration level is also obtained in the simulation, about 15 m away from the center of the tunnel   The frequencies are consistent with the findings by Kaewunruen that most vibration issues are between 2 and 150 Hz for the neighborhood near railway lines in tunnels, on embankments, in a cutting or on the level ground [14].In addition, there is a peak value within 5 Hz in the simulation, which is different from the measurement.The reason is that there is a fundamental frequency in the numerical model due to the limited element dimensions, and amplification effect generates as the excited frequency is close to the model's basic frequency.
are generally similar.The frequencies are consistent with the findings by Kaewunruen that most vibration issues are between 2 and 150 Hz for the neighborhood near railway lines in tunnels, on embankments, in a cutting or on the level ground [14].In addition, there is a peak value within 5 Hz in the simulation, which is different from the measurement.The reason is that there is a fundamental frequency in the numerical model due to the limited element dimensions, and amplification effect generates as the excited frequency is close to the model's basic frequency.5).In Figure 10, the calculated acceleration vibration levels are larger than that in the measurement as the 1/3 octave frequency is smaller than 5 Hz, which is mainly due to the amplification effect.On the other hand, the simulation vibration levels are consistent with the measurement in the higher frequencies, in which the trend of the 1/3 octave band frequency spectra by the observation and simulation results agree with the in-situ vibration measurements [7] and the observations on the waiting platform in Boston [45].Generally, the directivity differences in the measurement and simulation frequencies are small, demonstrating that the coupled finite element-infinite element boundary method can better capture the vibration characteristics in practical engineering, since the frequency of subway-induced vibration is high [43].5).In Figure 10, the calculated acceleration vibration levels are larger than that in the measurement as the 1/3 octave frequency is smaller than 5 Hz, which is mainly due to the amplification effect.On the other hand, the simulation vibration levels are consistent with the measurement in the higher frequencies, in which the trend of the 1/3 octave band frequency spectra by the observation and simulation results agree with the in-situ vibration measurements [7] and the observations on the waiting platform in Boston [45].Generally, the directivity differences in the measurement and simulation frequencies are small, demonstrating that the coupled finite element-infinite element boundary method can better capture the vibration characteristics in practical engineering, since the frequency of subway-induced vibration is high [43].

Simplified Prediction Formula of Ground Vibration Level
The ground vibration level is mainly influenced by the tunnel track-bed vibration and the tunnel depth.A simplified prediction formula of the ground vibration level of the track bed-tunnel liningsurrounding formation on this engineering certain soft site, based on the numerical simulation models with different tunnel depths and the track-bed accelerations, is presented as Equation ( 7).
( 0.0028 0.0033 0.942) where Z VL is the ground vertical acceleration level (dB), and ZS VL is the tunnel track-bed vertical acceleration level (dB).x is the horizontal distance away from the center of the tunnel (m), and h is the tunnel depth (m). Figure 11 plots the attenuation of peak vibration levels with respect to distance in this practical engineering.The predicted vibration levels are approximately similar to the measured results, with maximum error 3 dB.Therefore, the proposed prediction formula is acceptable and applicable for the prediction of vibration levels induced by subway in the practical engineering.

Simplified Prediction Formula of Ground Vibration Level
The ground vibration level is mainly influenced by the tunnel track-bed vibration and the tunnel depth.A simplified prediction formula of the ground vibration level of the track bed-tunnel lining-surrounding formation on this engineering certain soft site, based on the numerical simulation models with different tunnel depths and the track-bed accelerations, is presented as Equation (7).
where VL Z is the ground vertical acceleration level (dB), and VL ZS is the tunnel track-bed vertical acceleration level (dB).x is the horizontal distance away from the center of the tunnel (m), and h is the tunnel depth (m). Figure 11 plots the attenuation of peak vibration levels with respect to distance in this practical engineering.The predicted vibration levels are approximately similar to the measured results, with maximum error 3 dB.Therefore, the proposed prediction formula is acceptable and applicable for the prediction of vibration levels induced by subway in the practical engineering.

Parametric Study of Coupled Finite Element-Infinite Element Boundary Method in the Subway Induced Vibration
The elastic half-space simulation introduced in Section 2.2 is used in this section to discuss the efficiency of the presented finite element method coupled infinite element boundary and the parametric study is conducted to investigate the influence of model widths and depths on the numerical results.

The Effect of Model Depths on the Simulation Results
Table 2 shows the numerical cases with different model heights.The model can be expressed as LIN_H, where L denotes the width of the model, H is the height and IN means that the infinite element is adopted to the lateral side of the model.The base of the model is fixed in the simulation.The dimension of the element is 1 m × 1 m.Plain strain element CPE4 is used as the finite element in this model.Different model depths are analyzed to obtain a stable solution.Figure 12 is the sketch of the model with infinite element.Symmetrical constraints are imposed on the left boundary and the infinite element is applied to the right side of the model.

Parametric Study of Coupled Finite Element-Infinite Element Boundary Method in the Subway Induced Vibration
The elastic half-space simulation introduced in Section 2.2 is used in this section to discuss the efficiency of the presented finite element method coupled infinite element boundary and the parametric study is conducted to investigate the influence of model widths and depths on the numerical results.

The Effect of Model Depths on the Simulation Results
Table 2 shows the numerical cases with different model heights.The model can be expressed as LIN_H, where L denotes the width of the model, H is the height and IN means that the infinite element is adopted to the lateral side of the model.The base of the model is fixed in the simulation.The dimension of the element is 1 m × 1 m.Plain strain element CPE4 is used as the finite element in this model.Different model depths are analyzed to obtain a stable solution.Figure 12 is the sketch of the model with infinite element.Symmetrical constraints are imposed on the left boundary and the infinite element is applied to the right side of the model.

Parametric Study of Coupled Finite Element-Infinite Element Boundary Method in the Subway Induced Vibration
The elastic half-space simulation introduced in Section 2.2 is used in this section to discuss the efficiency of the presented finite element method coupled infinite element boundary and the parametric study is conducted to investigate the influence of model widths and depths on the numerical results.

The Effect of Model Depths on the Simulation Results
Table 2 shows the numerical cases with different model heights.The model can be expressed as LIN_H, where L denotes the width of the model, H is the height and IN means that the infinite element is adopted to the lateral side of the model.The base of the model is fixed in the simulation.The dimension of the element is 1 m × 1 m.Plain strain element CPE4 is used as the finite element in this model.Different model depths are analyzed to obtain a stable solution.Figure 12 is the sketch of the model with infinite element.Symmetrical constraints are imposed on the left boundary and the infinite element is applied to the right side of the model.The attenuation curves of the model applied to lateral infinite element boundary with different heights under the excitations of different vibration frequencies are shown in Figure 13.The numerical results of models with different heights are similar under the high frequency load, such as 40 Hz and 80 Hz, while the results are obviously different under the low frequency load.The numerical errors compared to the benchmark model 240IN_1280 mainly depend on the vertical frequencies, using an equation f 0 = c p /4H, with 3.8 Hz, 1.9 Hz, 0.95 Hz, 0.48 Hz and 0.24 Hz of the models in Cases 1-5, respectively.Therefore, the vertical frequencies of the models are small, resulting relative larger errors under low frequency load (close to the vertical fundamental frequency) than that under high frequency load.Vertical frequency, depending on the depth of the model with certain soil property, has great influence on the simulation results.Therefore, the depth of the model should be designed seriously to reduce the numerical errors due to the improper depth of the model.This agrees with the conclusion obtained by Andersen [46] that the accuracy of the simulation results increases with the distance from the tunnel; in other words, the depth of the model is significant to obtain more accurate estimates.
Figure 14 shows the amplitude-frequency curve of the models with different heights.It can be observed that: (1) There is an obvious peak in the curve when the model height is 80 m, while the amplitude-frequency curves are monotonic in other models.This is because the vertical frequency becomes smaller with the increase of model depth, and the vertical frequencies are less than 1 Hz to most models.However, the applied frequencies are equal to or greater than 1 Hz, larger than the vertical frequency.Therefore, there are no obvious peak values of most models in the amplitude-frequency curve.(2) The frequencies corresponding to the maximum error of the amplitude-frequency curves are mainly in the low frequency (1-6 Hz).Moreover, the error of amplitude-frequency curve is less than 5% when the model height is 640 m, indicating that numerical results are stable.The attenuation curves of the model applied to lateral infinite element boundary with different heights under the excitations of different vibration frequencies are shown in Figure 13.The numerical results of models with different heights are similar under the high frequency load, such as 40 Hz and 80 Hz, while the results are obviously different under the low frequency load.The numerical errors compared to the benchmark model 240IN_1280 mainly depend on the vertical frequencies, using an equation 0 / 4 p f c H  , with 3.8 Hz, 1.9 Hz, 0.95 Hz, 0.48 Hz and 0.24 Hz of the models in Cases 1-5, respectively.Therefore, the vertical frequencies of the models are small, resulting relative larger errors under low frequency load (close to the vertical fundamental frequency) than that under high frequency load.Vertical frequency, depending on the depth of the model with certain soil property, has great influence on the simulation results.Therefore, the depth of the model should be designed seriously to reduce the numerical errors due to the improper depth of the model.This agrees with the conclusion obtained by Andersen [46] that the accuracy of the simulation results increases with the distance from the tunnel; in other words, the depth of the model is significant to obtain more accurate estimates.
Figure 14 shows the amplitude-frequency curve of the models with different heights.It can be observed that: (1) There is an obvious peak in the curve when the model height is 80 m, while the amplitude-frequency curves are monotonic in other models.This is because the vertical frequency becomes smaller with the increase of model depth, and the vertical frequencies are less than 1 Hz to most models.However, the applied frequencies are equal to or greater than 1 Hz, larger than the vertical frequency.Therefore, there are no obvious peak values of most models in the amplitudefrequency curve.(2) The frequencies corresponding to the maximum error of the amplitudefrequency curves are mainly in the low frequency (1-6 Hz).Moreover, the error of amplitudefrequency curve is less than 5% when the model height is 640 m, indicating that numerical results are stable.

The Effect of Model Lateral Sides with Infinite Element Boundary on the Simulation Results
Table 3 shows the numerical cases of models with lateral infinite element.The applied boundary is similar to the models in Section 4.1.In this section, the model height remains unchanged to study the influence of model width on the numerical results with the infinite element applied to the lateral side.The numerical simulation shows that the solution is stable with the model width 1280 m, based on the discussion in Section 2.2.Therefore, the result of model 1280_80 is used as a benchmark to be compared with other models' results in this section.
The attenuation curves of the infinite element model with different widths under the excitation of vibration frequencies are shown in Figure 15.It can be seen that: (1) There is a certain requirement on the dimensions of the finite element area with the infinite element boundary applying to the lateral side, since smaller finite element area generates greater calculation error.While with the increasing of width in the finite element area, the calculation results tend to be stable gradually.(2) The simulation errors of finite element model 160_80 and model 160IN_80 are 60.78% and 13.88%, respectively, indicating that the maximum error is significantly reduced with the infinite boundary.The finite element model notes that infinite element boundary is not applied to the model.The relative error is the simulation results comparison with that in the benchmark model 1280_80.

The Effect of Model Lateral Sides with Infinite Element Boundary on the Simulation Results
Table 3 shows the numerical cases of models with lateral infinite element.The applied boundary is similar to the models in Section 4.1.In this section, the model height remains unchanged to study the influence of model width on the numerical results with the infinite element applied to the lateral side.The numerical simulation shows that the solution is stable with the model width 1280 m, based on the discussion in Section 2.2.Therefore, the result of model 1280_80 is used as a benchmark to be compared with other models' results in this section.
The attenuation curves of the infinite element model with different widths under the excitation of vibration frequencies are shown in Figure 15.It can be seen that: (1) There is a certain requirement on the dimensions of the finite element area with the infinite element boundary applying to the lateral side, since smaller finite element area generates greater calculation error.While with the increasing of width in the finite element area, the calculation results tend to be stable gradually.(2) The simulation errors of finite element model 160_80 and model 160IN_80 are 60.78% and 13.88%, respectively, indicating that the maximum error is significantly reduced with the infinite boundary.The finite element model notes that infinite element boundary is not applied to the model.The relative error is the simulation results comparison with that in the benchmark model 1280_80.Figure 16 shows the amplitude-frequency curve of the models with different widths and the infinite element boundary is applied to the lateral side in the simulation models.It can be observed that the errors of the infinite element model with 480 m width is less than 5%, and when the model width is 320 m, its corresponding error is 10%.The widths of the finite element model are 960 m and 640 m, respectively, to obtain the same accuracy.Therefore, the horizontal dimension of the finite element area can be reduced 50% effectively by means of the infinite element boundary, and the calculation accuracy remains unchanged.Moreover, the computational time can be saved due to the decreasing number of element units, demonstrating the higher performance of the combined finite element-infinite element boundary method.Figure 16 shows the amplitude-frequency curve of the models with different widths and the infinite element boundary is applied to the lateral side in the simulation models.It can be observed that the errors of the infinite element model with 480 m width is less than 5%, and when the model width is 320 m, its corresponding error is 10%.The widths of the finite element model are 960 m and 640 m, respectively, to obtain the same accuracy.Therefore, the horizontal dimension of the finite element area can be reduced 50% effectively by means of the infinite element boundary, and the calculation accuracy remains unchanged.Moreover, the computational time can be saved due to the decreasing number of element units, demonstrating the higher performance of the combined finite element-infinite element boundary method.Figure 16 shows the amplitude-frequency curve of the models with different widths and the infinite element boundary is applied to the lateral side in the simulation models.It can be observed that the errors of the infinite element model with 480 m width is less than 5%, and when the model width is 320 m, its corresponding error is 10%.The widths of the finite element model are 960 m and 640 m, respectively, to obtain the same accuracy.Therefore, the horizontal dimension of the finite element area can be reduced 50% effectively by means of the infinite element boundary, and the calculation accuracy remains unchanged.Moreover, the computational time can be saved due to the decreasing number of element units, demonstrating the higher performance of the combined finite element-infinite element boundary method.

The Effect of Model Base with Infinite Element Boundary on the Simulation Results
Table 4 shows the numerical cases with different model depths.The model can be expressed as L IN_H IN.The infinite element is applied to the right side and the base of the model.In these cases, the model width remains unchanged to study the influence of model depth on the numerical results with the infinite element.The numerical simulation shows that the solution is stable when the model's height is 1280 m based on the discussion in Section 4.1.Therefore, the result of model 240IN_1280 is used as a benchmark to be compared with other model results in this section.The attenuation curves of the infinite element models with different heights under the excitation of vibration frequencies are shown in Figure 17, from which it can be concluded that there is a relatively large difference between the calculated curve and the stable solution under the low frequency excitation when the model depth is small.However, the results are close to the stable solution under the high frequency excitation, such as 40 Hz.Therefore, larger height and higher frequency load generally generate smaller simulation errors.

The Effect of Model Base with Infinite Element Boundary on the Simulation Results
Table 4 shows the numerical cases with different model depths.The model can be expressed as L IN_H IN.The infinite element is applied to the right side and the base of the model.In these cases, the model width remains unchanged to study the influence of model depth on the numerical results with the infinite element.The numerical simulation shows that the solution is stable when the model's height is 1280 m based on the discussion in Section 4.1.Therefore, the result of model 240IN_1280 is used as a benchmark to be compared with other model results in this section.The attenuation curves of the infinite element models with different heights under the excitation of vibration frequencies are shown in Figure 17, from which it can be concluded that there is a relatively large difference between the calculated curve and the stable solution under the low frequency excitation when the model depth is small.However, the results are close to the stable solution under the high frequency excitation, such as 40 Hz.Therefore, larger height and higher frequency load generally generate smaller simulation errors.

The Effect of Model Base with Infinite Element Boundary on the Simulation Results
Table 4 shows the numerical cases with different model depths.The model can be expressed as L IN_H IN.The infinite element is applied to the right side and the base of the model.In these cases, the model width remains unchanged to study the influence of model depth on the numerical results with the infinite element.The numerical simulation shows that the solution is stable when the model's height is 1280 m based on the discussion in Section 4.1.Therefore, the result of model 240IN_1280 is used as a benchmark to be compared with other model results in this section.The attenuation curves of the infinite element models with different heights under the excitation of vibration frequencies are shown in Figure 17, from which it can be concluded that there is a relatively large difference between the calculated curve and the stable solution under the low frequency excitation when the model depth is small.However, the results are close to the stable solution under the high frequency excitation, such as 40 Hz.Therefore, larger height and higher frequency load generally generate smaller simulation errors.Figure 18 shows the amplitude-frequency curve of the models with different heights as infinite element boundary is applied to the lateral side and base of the model.It can be observed that the results of amplitude-frequency in the infinite element models are much closer to the stable solution compared with that in Figure 14 in Section 4.1.The reason is that the fundamental frequency of the infinite models decreases, close to that of benchmark 240IN_1280.In addition, from the energy point of view, the reflection energy becomes weaker and close to the half-space infinite body after applying the infinite-element boundary.The near-source induced waves usually contain high energy, which may cause greater vibration [44].
Moreover, the model 240IN_640 in Section 4.1 can be replaced by the model 240IN_240IN in this section to achieve the same calculation accuracy, demonstrating that the height of finite element model can be reduced 62.5% after applying the infinite element boundary to the model base.Meanwhile, computational time is saved effectively, revealing the high performance of the coupled finite element-infinite element boundary method in decreasing the models' dimensions once more.Figure 18 shows the amplitude-frequency curve of the models with different heights as infinite element boundary is applied to the lateral side and base of the model.It can be observed that the results of amplitude-frequency in the infinite element models are much closer to the stable solution compared with that in Figure 14 in Section 4.1.The reason is that the fundamental frequency of the infinite models decreases, close to that of benchmark 240IN_1280.In addition, from the energy point of view, the reflection energy becomes weaker and close to the half-space infinite body after applying the infinite-element boundary.The near-source induced waves usually contain high energy, which may cause greater vibration [44].
Moreover, the model 240IN_640 in Section 4.1 can be replaced by the model 240IN_240IN in this section to achieve the same calculation accuracy, demonstrating that the height of finite element model can be reduced 62.5% after applying the infinite element boundary to the model base.Meanwhile, computational time is saved effectively, revealing the high performance of the coupled finite element-infinite element boundary method in decreasing the models' dimensions once more.Figure 18 shows the amplitude-frequency curve of the models with different heights as infinite element boundary is applied to the lateral side and base of the model.It can be observed that the results of amplitude-frequency in the infinite element models are much closer to the stable solution compared with that in Figure 14 in Section 4.1.The reason is that the fundamental frequency of the infinite models decreases, close to that of benchmark 240IN_1280.In addition, from the energy point of view, the reflection energy becomes weaker and close to the half-space infinite body after applying the infinite-element boundary.The near-source induced waves usually contain high energy, which may cause greater vibration [44].
Moreover, the model 240IN_640 in Section 4.1 can be replaced by the model 240IN_240IN in this section to achieve the same calculation accuracy, demonstrating that the height of finite element model can be reduced 62.5% after applying the infinite element boundary to the model base.Meanwhile, computational time is saved effectively, revealing the high performance of the coupled finite element-infinite element boundary method in decreasing the models' dimensions once more.

Conclusions
In this study, a coupled finite element-infinite element boundary model based on the basic principles of the infinite element method was performed by the harmonic analysis based on the commercial software ABAQUS to address the issue of subway induced vibration, a major environmental concern in urban areas.In addition, the finite element method in the harmonic analysis was verified by comparing it with the thin layer method.Moreover, an interface program was developed to automatically read the simulation result files (dat file) in the ABAQUS harmonic analysis, and then the data were put into the frequency domain calculation on the MATLAB platform.Besides, a measurement of a practical engineering was performed on a track bed-tunnel lining-surrounding formation located on Line 2, Shanghai Metro, and the corresponding simulation model was established based on the introduced method.In this paper, the numerical results of vibration level and spectrum characteristic are discussed.Moreover, a simplified prediction formula of ground vibration level is proposed based on the measurement results.Finally, parametric study was conducted in an elastic half-space simulation to investigate the influence of model widths and depths with infinite element boundary on the numerical results, demonstrating the effectiveness and high performance of the coupled finite element-infinite element boundary model in the harmonic analysis based on ABAQUS.The following are the conclusions drawn from this study, which may provide insight towards the rational models of subway induced vibration.

1.
The coupled simulation method is applicable for the practical engineering projects of subway induced vibration on the soft site in Shanghai, and the proposed formula for the vibration levels is reasonable for the prediction in the subway induced vibration.

2.
There exists a vibration amplifying zone a certain distance from the tunnel under high frequency loads due to wave propagation and reflection, and the amplification has relationship with the spectrum of the applied excitation accelerations.

3.
The finite element model dimensions can be reduced 50% and 66.67% effectively in the widths and depth, respectively, with the coupled method, applying the infinite element boundary to the lateral and base sides of the finite element model.

4.
The model dimensions have great influence on the numerical results as the excited frequency is close to the model's vertical fundamental frequency, and the depth of the model should be designed seriously to reduce the numerical errors.

Figure 2 .
Figure 2. The attenuation curves of TLM results and ABAQUS results with two different element dimensions.Note: 1 and 0.25 mean the element dimensions are 1 m and 0.25 m, respectively.TLM is the thin layer method.

Figure 3 .
Figure 3. (a) The simulation model; and (b) the applied displacement load time-history at Point A.

Figure 4 .
Figure 4. (a) Vertical displacement time-history; and (b) local displacement time-history of Point B.

FieldFigure 3 .
Figure 3. (a) The simulation model; and (b) the applied displacement load time-history at Point A.

Figure 3 .
Figure 3. (a) The simulation model; and (b) the applied displacement load time-history at Point A.

Figure 4 .
Figure 4. (a) Vertical displacement time-history; and (b) local displacement time-history of Point B.
The acceleration time-histories of the central tunnel and the typical ground points were measured and acquired in this test.The track bed-tunnel lining-surrounding formation model and the arrangement of measuring points in lateral direction, which were 0 m, 5 m, 11 m, 15 m, 25 m and 40 m from the central tunnel, are shown in Figure 5.The outer radius of the tunnel is 3.1 m and inner radius is 2.75 m.The tunnel depth, the distance from the top of the tunnel to the surface of the tunnel, is 13.8 m.The elastic modulus, the density and

Figure 4 .
Figure 4. (a) Vertical displacement time-history; and (b) local displacement time-history of Point B.

Figure 5 .Figure 5 .
Figure 5.The schematic diagram of the model and the arrangement of measuring points.

Figure 5 .Figure 6 .
Figure 5.The schematic diagram of the model and the arrangement of measuring points.

Figure 6 .
Figure 6.(a) The acceleration time-history of the track-bed; (b) the corresponding Fourier spectrum of the track-bed; (c) 1/3 octave frequency; and (d) acceleration levels of track-bed.

20 Figure 6 .
Figure 6.(a) The acceleration time-history of the track-bed; (b) the corresponding Fourier spectrum of the track-bed; (c) 1/3 octave frequency; and (d) acceleration levels of track-bed.

Figure 7 .
Figure 7.The meshed model of track bed-tunnel lining-surrounding formation.

Figure 7 .
Figure 7.The meshed model of track bed-tunnel lining-surrounding formation.

Figure 8 .
Figure 8.(a) Vibration levels with respect to distance under the 31 measured accelerations excitations; (b) vibration levels in Cases 11 and 29; and (c) the mean vibration levels of measurement and simulation.Note: In (c), M and S mean the measurement results and the simulation results.

Figure 9 2 )Figure 8 .
Figure9plots the simulation and measurement frequency spectra 0 m, 5 m and 15 m from the center tunnel.The calculated and measured frequency spectra, with dominant frequency 40-90 Hz, are generally similar.The frequencies are consistent with the findings by Kaewunruen that most vibration issues are between 2 and 150 Hz for the neighborhood near railway lines in tunnels, on embankments, in a cutting or on the level ground[14].In addition, there is a peak value within 5 Hz in the simulation, which is different from the measurement.The reason is that there is a fundamental frequency in the numerical model due to the limited element dimensions, and amplification effect generates as the excited frequency is close to the model's basic frequency.

Figure 9
Figure9plots the simulation and measurement frequency spectra 0 m, 5 m and 15 m from the center tunnel.The calculated and measured frequency spectra, with dominant frequency 40-90 Hz, are generally similar.The frequencies are consistent with the findings by Kaewunruen that most vibration issues are between 2 and 150 Hz for the neighborhood near railway lines in tunnels, on embankments, in a cutting or on the level ground[14].In addition, there is a peak value within 5 Hz in the simulation, which is different from the measurement.The reason is that there is a fundamental frequency in the numerical model due to the limited element dimensions, and amplification effect generates as the excited frequency is close to the model's basic frequency.

Figure 10
Figure 10 displays the compassion 1/3 octave band frequency spectra of simulation with the letter S drawn by the blue lines and measurement with the letter M to show the trend of the vertical response at the distance of 0 m, 5 m, 11 m and 15 m to the central tunnel (as shown in Figure5).In Figure10, the calculated acceleration vibration levels are larger than that in the measurement as the 1/3 octave frequency is smaller than 5 Hz, which is mainly due to the amplification effect.On the other hand, the simulation vibration levels are consistent with the measurement in the higher frequencies, in which the trend of the 1/3 octave band frequency spectra by the observation and simulation results agree with the in-situ vibration measurements[7] and the observations on the waiting platform in Boston[45].Generally, the directivity differences in the measurement and simulation frequencies are small, demonstrating that the coupled finite element-infinite element boundary method can better capture the vibration characteristics in practical engineering, since the frequency of subway-induced vibration is high[43].

Figure 10 displays
Figure 10 displays the compassion 1/3 octave band frequency spectra of simulation with the letter S drawn by the blue lines and measurement with the letter M to show the trend of the vertical response at the distance of 0 m, 5 m, 11 m and 15 m to the central tunnel (as shown in Figure5).In Figure10, the calculated acceleration vibration levels are larger than that in the measurement as the 1/3 octave frequency is smaller than 5 Hz, which is mainly due to the amplification effect.On the other hand, the simulation vibration levels are consistent with the measurement in the higher frequencies, in which the trend of the 1/3 octave band frequency spectra by the observation and simulation results agree with the in-situ vibration measurements[7] and the observations on the waiting platform in Boston[45].Generally, the directivity differences in the measurement and simulation frequencies are small, demonstrating that the coupled finite element-infinite element boundary method can better

Figure 11 .
Figure 11.The comparison of vibration levels of the engineering project by the measurement and presented prediction formula.

Figure 11 .
Figure 11.The comparison of vibration levels of the engineering project by the measurement and presented prediction formula.

Figure 11 .
Figure 11.The comparison of vibration levels of the engineering project by the measurement and presented prediction formula.

Figure 12 .
Figure 12.The diagram of finite element method coupled infinite element boundary.

Figure 14 .
Figure 14.The amplitude-frequency curve with different model depths at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Figure 14 .
Figure 14.The amplitude-frequency curve with different model depths at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Figure 16 .
Figure 16.The amplitude-frequency curves of models with lateral infinite element boundary at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Figure 16 .
Figure 16.The amplitude-frequency curves of models with lateral infinite element boundary at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Figure 18 .
Figure 18.The amplitude-frequency curves of models with base infinite element boundary at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Figure 18 .
Figure 18.The amplitude-frequency curves of models with base infinite element boundary at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Figure 18 .
Figure 18.The amplitude-frequency curves of models with base infinite element boundary at the ground evaluation point of: (a) 2 m; (b) 10 m; (c) 16 m; and (d) 30 m.

Table 1 .
Soil parameters of different soil layers.
in Section 2.3.
Poisson's ratio of concrete are 34,500 MPa, 2500 kg/m 3 and 0.25, respectively.The typical clay in Shanghai, not homogeneous in the depth direction, was employed in this field test and numerical simulation.The soil parameters of different soil layers are shown in Table1in Section 2.3.

Table 2 .
Numerical cases with different model depths.

Table 2 .
Numerical cases with different model depths.

Table 2 .
Numerical cases with different model depths.

Table 3 .
Numerical cases of different model widths with lateral infinite element.

Table 3 .
Numerical cases of different model widths with lateral infinite element.

Table 4 .
Numerical cases of different model depths with base infinite element.

Table 4 .
Numerical cases of different model depths with base infinite element.

Table 4 .
Numerical cases of different model depths with base infinite element.