Abstract
With the fast development of rail transit, the environmental vibration problems caused by subways have received increasing attention. A 3D finite element model was built in this study to investigate the ground vibrations induced by the moving load operating in the parallel twin tunnels. Compared to the model consisting of a single tunnel that was commonly adopted in the past studies, a pair of tunnels is considered and the surrounding medium of the tunnels is taken as a saturated porous medium. The governing equations of the 3D finite element method modeling of the saturated poroelastic soil have been derived according to Biot’s theory. Computed results showed that the dynamic response of the twin-tunnel model is greater than that of the single tunnel model. And the spacing between two tunnels, tunnel buried depth, and load moving speed are the essential parameters to determine the dynamic response of the tunnel and soil.
1. Introduction
Recently, concerns about the environmental vibrations induced by the rail transit have increased substantially. When the subways run close to the existing infrastructures, vibrations are transmitted to the infrastructure through the ground, which can cause annoyance to inhabitants or result in malfunction interruption to sensitive equipment. Due to the high groundwater table in southeastern area of China, e.g., Shanghai, Zhejiang, Jiangsu, and Guangdong, subways in these areas are running in the saturated soil. Besides it is a common practice to construct underground railway lines in pairs. The vibrations caused by twin tunnels cannot be simplified to be the sum of those caused by two single tunnels. Therefore, an investigation on ground vibrations induced by subways in the context of twin tunnels embedding in a saturated ground is necessary and desirable.
The vibrations caused by subways have been investigated by many scholars using both analytical and computational methods. Metrikine et al. [1] presented an analytical method to simplify the metro to Euler beams embedded in viscoelastic medium. Their study focused on the effect of load velocity on the structural response. Forrest and Hunt [2, 3] built a three-dimensional analytical model for a circular subway tunnel buried deep underground in order to obtain its dynamic responses. Liu et al. [4, 5] studied the transient response of partially sealed spherical cavity embedded in viscoelastic saturated soil. The partial permeability of the boundary and the influence of the relative stiffness between lining and soil on the transient response were investigated. Yuan et al. [6] analyzed the influences of the load velocity and the oscillating frequencies on the structural displacement of the track, the displacement, and pore pressure responses of the ground. Considering the train load excitation as a random process, Hunt et al. [7] studied the ground vibrations due to the subway trains using a track-building model of infinite length.
Compared with the analytical method, the computational method, e.g., finite element method (FEM), stands out for its easy handling of the material heterogeneity and complex/irregular geometry. Balendra et al. [8] investigated the coupled vibrations between the subway, the ground, and the building using a 2D plane-strain finite element model. A viscous absorbing boundary condition was applied to the truncation boundary to suppress unwanted reflections. A 2D finite element model has been established by Wang et al. [9] including the subway tunnel, the ground, and the adjacent building. The impact of the tunnel depths on the vibrations of different floors and positions in the buildings was investigated. Yang et al. [10] studied the dynamic interaction between the ground and tunnel system using the finite element method. At the truncation boundary, the infinite element was applied to damp out the outgoing waves. Similarly, Zhang and Pan [11] established a coupled finite and infinite element model to study the vibration of metro tunnel structure and surrounding stratum. A coupled finite element (FE)–boundary element (BE) approach was proposed by Jones et al. [12] to study the wave transmissions and the ground vibrations induced by the subways.
However, the wave field generated by the moving load is three-dimensional in nature, which restricts the applicability of 2D computational models. Gardien and Stuit et al. [13] analyzed the ground vibrations generated by the subway using the finite element code LS-DYNA. Wolf et al. [14] built a 3D prediction model for low-frequency ground vibration induced by the subway using the finite difference method, which is calibrated using the empirical relationship established on the measured field data. G. Degrande et al. [15] proposed a coupled finite element (FE)–boundary element (BE) formulation to study the dynamic interaction between the tunnel and the ground, in which the tunnel was modeled by FEM and the ground by BEM. The dynamic responses of a shallow cut-and-cover masonry tunnel were compared to those of a deep bored one. A 3D periodic FE-BE model was established by Liu et al. [16] to obtain the train-induced vibrations for both the tunnel and the free field. The floating slab track and the general slab track are compared for their effectiveness in vibration reduction.
As mentioned above, the current research focuses on the environmental vibrations problems generated by subways running in a single tunnel. The effect of the neighboring tunnel is not considered although it is a common practice that underground railway lines are constructed in pairs. The vibration caused by two metro tunnels is different to the sum of two single tunnels since the neighboring tunnel can impede and screen the vibration waves generated by the operating tunnel. To the best of the authors’ knowledge, Kuo et al. [17] presented the investigation on the dynamics of a twin tunnel embedded in elastic full-space, adopting a coupled finite element–boundary element approach. Parametric studies were carried out for the load excitation frequency, the tunnel orientations, and the tunnel geometry, which demonstrated that the dynamic interactions between the two tunnels were highly significant. However, the full space model used by Kuo et al. is not applicable when the ground vibrations are of concern, since the ground is semi-infinite with a free surface. Moreover, the possible influences of the load speed, the tunnel burying depth, and the tunnel spacing on the dynamics responses of the twin-tunnel and the ground remain untouched in the literature.
The ground in the above-referred papers using numerical approach is generally modeled as single-phase elastic medium in which the ground water is not considered. Existing research by Yuan et al. [18] showed that the dynamic responses of the subway tunnel in the saturated ground are different from those in the elastic ground. Therefore, it is necessary to model the ground as saturated poroelastic medium such that the ground water can be taken into account.
This paper focuses on the coupled vibrations of the liner, the twin-parallel tunnels, and the saturated ground using the finite element method. A moving point load is applied to the invert to represent the subway excitation. The dynamics of the surrounding saturated soil are governed by Biot’s poroelastic theory. To suppress spurious reflections at the truncation boundaries, an artificial boundary condition named the multi-transmitting formula (MTF) that is proposed by the Liao and Wong [19] and extended by Shi et al. [20] is applied to transmit the outgoing waves in the saturated soil medium. The effectiveness and stability of MTF have been demonstrated in the paper by Shi et al. [20]. Using the established 3D finite element model, parametric studies have been performed for investigating the effects of the load velocity, the burying depth of tunnel, and the tunnel spacing on the coupled vibrations of the liner-tunnel-ground system. Because the effect of the permeability coefficient on the dynamic response of soil has been fully discussed by many researchers and we have not found any new conclusions using the presented model, therefore, the influence of the permeability coefficient on vibration response has not been discussed in this paper.
2. Governing Equations for Saturated Ground
According to Biot’s theory [21, 22], the governing equations for the dynamics of the saturated poroelastic medium read as follows:where denotes the gravity force acting on the soil skeleton; and denote the density of the soil and the density of the pore fluid, respectively; , is the density of the solid skeleton and is the porosity; and denote the displacement of soil skeleton and the infiltration displacement of the pore fluid relative to the soil skeleton, respectively; “” and “” denote the first and second derivative of time, respectively; ; , is the acceleration of gravity, is the Darcy permeability; and denote the total and effective stress tensors, respectively; is the Kronecker delta; accounts for the pore pressure; and are Biot parameters; and subscripts i and j denote the tensor operation and the summation convention is applied.
After manipulation, the pore pressure can be eliminated from the governing equations, leaving and the field variables:
There are two reasons for keeping and as the primary unknowns: the first reason is that the formulation is a complete description of Biot’s theory. For the formulation, the relative pore fluid acceleration needs to be neglected [23], which would bring in certain inaccuracy especially for a high-frequency loading situation [23, 24]; the second reason is that the formulation is consistent with the proposed artificial boundary condition, multi-transmitting formula (MTF, given in Section 3), since it is vector-based and manipulates the displacement vectors and directly, while the pore pressure p is a scalar and thus not a suitable variable for applying the MTF.
The same interpolation function is used for both the soil-skeleton displacement and the relative pore-fluid displacement which are given aswhere is the coordinate; represents time; is the shape function at node ; is the number of nodes of each element; is the unit matrix of 2 × 2 or 3 × 3; and and are the nodal displacement vectors.
Following the standard Galerkin procedure, the governing equations of the finite element formulation representing the saturated poroelastic soil medium can be derived as
The expressions of the elements in the mass matrix, the damping matrix, the stiffness matrix, and the force vector are given in the appendix.
3. The Multi-Transmitting Formula
Due to the semi-infinite nature of the ground, the ground vibration analysis using the finite element method requires artificial boundary conditions to make the computational domain finite. Many artificial boundary conditions have been proposed for absorbing/damping the outgoing waves for both the elastic medium and the saturated poroelastic medium. For a detailed review, one is referred to the paper by Shi et al. [20]. Among the existing artificial boundary conditions, the multi-transmitting formula, that is proposed by Liao and Wong [19] for the elastic medium and then extended by Shi et al [20] for the saturated poroelastic medium, stands out for being local in both time and space and for the easy implementation in the finite element analysis both for 2D and 3D grids.
The MTF extrapolates displacement on the artificial boundary at time as a linear combination of the displacements at previous time steps along a straight line normal to the boundary, which is given as
where denotes the soil skeleton displacements of the grid point at and at time . Similarly, denotes the relative pore-fluid displacement; represents the apparent velocities of waves propagating along the normal direction of artificial boundary; is the binomial coefficient; and N is the order of the boundary condition.
Since the computational points () in (9) do not generally coincide with the grid points, a quadratic interpolation scheme is employed to correlate the displacements at the computational points to those at the grid points aswhere the operator A is the backward space shift of :in which, the interpolation coefficients are given as
The operators and are the forward time and space shifts, respectively:
Correspondingly, the operators and represent the backward time and space shifts, respectively.
According to Shi et al. [20], a second-order MTF () and the artificial wave velocity ( is the shear wave velocity of the saturated medium) are efficient and accurate enough for configuring the MTF. Moreover, the numerical stability of the MTF can be guaranteed by (1) meeting the Courant-Friedrichs-Lewy (CFL) condition and (2) perturbing the interpolation coefficient by adding a small negative value ().
4. Verification of the Present Model
To verify the accuracy of the proposed finite element formulation and the effectiveness of the MTF boundary condition, the transient response due to an edge pressure applied radially at a cavity surface is computed and compared to the analytical solution given by T. Senjuntichai [25], which is developed for a cavity in an infinite saturated soil medium under the plane-strain condition.
Due to the symmetry, a 1/4 model is established for the cavity and the surrounding saturated soil medium. Symmetric boundary conditions are imposed for the left and bottom boundaries, while MTF boundary conditions are applied to the top and right boundaries to account for the infinite extension of the soil medium. The mesh in the x-y plane is shown in Figure 1(a), while the movement in the z direction is fixed to model the plain strain condition. The applied pressure is a triangular pulse, as given in Figure 1(b), in which is the normalized time , is the radius of tunnel, denotes the edge pressure, and is its peak value.
(a)
(b)
The soil properties and the load parameters are the same as those given by Senjuntichai [25] and summarized in Table 1. Figure 2 shows the comparison on the radial displacement at the cavity surface, where the vertical axis is the normalized radial displacements and the horizontal axis is the normalized time . From Figure 2, certain reflections can be observed for the MTF condition when the reflected waves reached the observation point at the first time (i.e., ). However, the following impinging of the reflected waves is transmitted effectively by MTF rending a very close comparison to the analytical solution for . Also it can be seen from Figure 2 that significant reflections can happen during the entire time if the fixed boundary condition is used.
5. Numerical Results
In this section, coupled vibrations of the liner, the twin tunnels, and the saturated ground are investigated using the developed finite element along with the developed MTF boundary condition. Parametric studies have been performed for the burying depth of tunnel, the tunnel spacing, the load velocity, and the soil permeability.
5.1. Model Description
The computation model consisting of two identical circular tunnels that are parallel to each other and embedded in a three-dimensional saturated ground is shown in Figure 3. The model size is (width) × (height) × (length). The center-to-center spacing between the tunnels is represented by , and the burying depth measured from the ground surface to the tunnel center is denoted by , as shown in Figure 3(b). denotes the outer radius of the liner and the thickness of liner is represented by . The side and cross-sectional views with dimensions are presented in Figures 3(a) and 3(b), respectively. A nonoscillating point load of magnitude that moves along the positive axle acts vertically to the invert of the left tunnel liner.
(a) Side view
(b) Cross section view
Four observation points are chosen to record the responses caused by the moving load: points A and B are located at the liner bottom of the tunnel with the load applied and the other tunnel, respectively; points C and D are located on the ground surface right between the twin tunnels and above the center of the tunnel loaded, respectively. All observation points are located at the middle of the model along the tunnel axial direction (i.e., ).
The bottom surface of the model is fixed to consider the underlying hard stratum. The top surface is free and set as permeable. The MTF boundary condition is applied to the remaining 4 side surfaces to account for the infinite extension of the saturated ground. For a reliable modeling of wave propagations, the element sizes of the ground should meet the requirement , where is the shear wave length and is the highest frequency of the shear wave [26]. The implicit Newmark method is employed for the time marching of the dynamic analysis with the Newmark parameter being set to 0.6 to remove disturbances caused by sudden application of the point load. Although the Newmark integration scheme is unconditionally stable, the time step should meet the CFL condition as mentioned in Section 3, since the MTF is essentially a forward Euler scheme.
Parameters of the liner, the saturated ground, and the moving load are summarized in Table 2. In the subsequent analysis, the dynamic displacement and velocity are presented in decibels (dB); i.e., , where denotes an arbitrary variable and is the reference value. The reference values for displacement and velocity are 2.0 × 10−5 m and 2.0 × 10−5 m/s, respectively.
5.2. Models Comparisons
In this section, a single tunnel model and a model with both tunnels loaded are compared with the present twin-tunnel model. The cross sections of the two models are presented in Figures 4(a) and 4(b), respectively. The results of single phased model are also compared with those of the saturated poroelastic model in this section.
(a) Single tunnel model
(b) Model with both tunnels loaded
To investigate the influence of interactions between the two tunnels on the vibration-prediction results, the comparison of dynamic responses between single tunnel model and the present model is shown in Figures 5(a) and 5(b). The burying depth and the spacing of the tunnels are, respectively, fixed at = 4 and = 4, and the load velocity is . In Figure 5(a), the tunnel spacing for the single tunnel model denotes two times the distance between tunnel center and the model center. The dynamic responses between two models on point A are very close while the results of the present model are greater than those of the single tunnel model on the ground observation points (point C and point D). The displacements at point C for the single tunnel model decrease faster than the present model between which is 4 to 5. That is because the dynamic response at point C in the present model is the sum of the waves generated from the tunnel with the load applied and reflected wave from the other tunnel lining. However, the dynamic responses of the single-tunnel model only come from the tunnel with the load applied. As the tunnel spacing increases, the differences of the displacement at point C and D for both models decrease. In Figure 5(b), it is observed that the maximum displacements of single tunnel model are all smaller than those of the present model when the burying depth changes from 3 to 5. The decrease rates for both models are almost the same.
(a) Maximum displacement at the observation points for different tunnel spacing
(b) Maximum displacement at the observation points for different burying depth
The comparison between the single tunnel model and the present model is shown in Figure 6 under different load velocities. It is observed from Figure 6 that the dynamic responses of the twin-tunnel model are larger than those of the single tunnel model which is inconsistent with the results shown in Figure 5. And the difference of the two models at point C is greater than that at point D. The differences reach peak value when . Since the dynamic response at point C is jointly affected by two tunnels while the dynamic response at point D is dominated by the tunnel loaded, therefore it is essential to consider the twin-tunnel model to obtain more accurate prediction results.
Practically, both metro tunnels are always loaded at the same time due to the high operation frequency of the subways. The dynamic responses of the model with both tunnel inverts loaded are studied in this section.
In contrast to the result of the present model, the dynamic responses of point C and point D due to two moving point loads are larger in terms of the maximum displacement. In Figure 7, the maximum displacements of the two models at point A are almost the same. The dynamic response of the two models at point C decreases a little with increasing tunnel spacing. However, as tunnel spacing increases, the maximum displacement of the two models at point A remains the same and the results of two models become closer at point D. This is because as the tunnel spacing increases, the effect of the other tunnel on point D decreases gradually. And with the increasing of the burying depth, the dynamic responses of points C and D have obvious downtrend. Besides, the displacement at the point C is larger than that at the point D in the model with two loads while the situation of the present model is the opposite. That is because the load in the other tunnel has a larger contribution to the dynamic response at point C than at point D.
(a) Maximum displacement at the observation points for different tunnel spacing
(b) Maximum displacement at the observation points for different burying depth
The comparison of the displacements and the vibration velocity between the present model and the model with both tunnels loaded is shown in Figures 8 and 9 for different load velocities. For points C and D, the dynamic responses of the model with both tunnels loaded are obviously larger than the model with only one tunnel loaded. However, the prediction results of the two models at point A are very close. And the maximum displacement at point C is smaller than that at point D in the model with only one tunnel loaded while the situation of the two-load model is the opposite which is explained in the previous section. The difference between the maximum displacements of the two models is almost unchanged at different load velocity. However, the difference in the vibration velocity results of the two models at point D decreases as the velocity increases and reaches a minimum at = 0.9. From Figure 9, it can also be observed that the maximum vibration velocities at points C and D are almost the same in the present model with only one tunnel loaded, while the value at point C in the two-load model is significantly larger than that at point D. Because the load in the other tunnel is closer to point C than point D, it has a greater impact on the dynamic response of point C.
To illustrate the necessity in introducing the saturated ground model, the dynamic responses of the twin tunnel embedded in a single-phase ground are compared to those of the saturated ground model and presented in Figure 10. The single-phase ground model is obtained by reducing the saturated ground model after fixing the pore-fluid relative displacement and setting Biot’s parameters to zeros. It is observed that the dynamic responses of the elastic ground are larger than those in the saturated ground. The discrepancy between the two ground models widens when the burying depth increases. And as the load velocity increases, the difference of dynamic response between the two models at point D also becomes larger. From the comparison it can be concluded that the ground water has an obvious influence on the vibration responses of the tunnel and the ground. Their vibration levels would be overestimated if the ground water is not considered.
(a) Maximum vibration velocity at observation points for different burying depth
(b) Maximum displacement at observation points for different load velocity
(c) Maximum vibration velocity at observation points for different load velocity
5.3. The Influences of Load Velocity
The influences of load velocity on the dynamic responses of the liner-tunnel-ground system are investigated in this section. Six different velocities = 0.1, 0.3, 0.5, 0.7, 0.9, and 1.0 are considered along with the tunnel spacing = 4 and the burying depth = 4.
It is observed from Figure 11 that the displacements at points A and B change slightly with the increasing load velocity, while the displacements at points C and D reach a peak at around . Since C and D are on the ground surface, their responses will be amplified when the load velocity approaches the Rayleigh wave velocity of the ground, which is around . The same results were reported in previous studies on the critical velocity for a train running on the ground surface [27, 28]. However, points A and B are set on the concrete liner, whose shear wave velocity is much higher than that of the soil. Thus the load velocity is very low compared to the shear wave velocity of the concrete which makes the displacement of points A and B nonsensitive to the change of load velocity.
5.4. The Influence of Tunnel Spacing
In this section, the influence of the tunnel spacing on the dynamic responses is investigated by considering 6 different twin tunnel distances, i.e., = 3, 4, 5, 6, 7, 8. The load velocity is . The maximum vertical displacements at the observation points when = 4 are compared in Figure 12.
The displacement at point A is significantly larger than those at other points, since point A is traversed directly by the point load. When the tunnel spacing increases, the displacements at points B and C decrease gradually. However, the displacements at points A and D remain almost the same. This observation is within expectation because the positions of points A and D remain the same while points B and C are shifted away from the load when the spacing increases. The displacements at points B and C vary greatly when the tunnel spacing increases from 4 to 6.
5.5. The Influence of Burying Depth
In this section, parametric study is conducted for 5 different burying depths of the twin tunnel, i.e., = 2, 3, 4, 5, 6. And the load velocity is set to .
The variation of the maximum vertical displacements at the observation points with respect to the burying depth is presented in Figure 13. In this case, the tunnel spacing is fixed as = 4. It is observed that the displacements at points A, C, and D decrease when increases. For the displacement at point B, it increases slightly when < 4; however, it drops substantially when increases further to 6. This phenomenon is explained by the fact that when the tunnel is shallowly buried, the waves generated by the source tunnel (i.e., the tunnel where the moving load is applied) impinge on the ground surface at a large incident angle (resembling the glancing incident); correspondingly, the reflected waves leave the surface at a large angle, and thus the surface reflection region at nonsourced tunnels is small resulting in a fraction of surface reflections contributing to the dynamic response at point B. However, the surface reflection region will expand when h/r increases gradually; but when the burying depth is large enough (i.e., > 5), the material and geometrical damping are dominant, which results in the decreasing of the displacement at point B.
6. Conclusions
In this paper, a 3D twin-tunnel model is presented to evaluate the dynamic responses of saturated half-space induced by subways. The effects of spacing between two tunnels, tunnel buried depth, load moving speed, and soil permeability on the vibration response are investigated. Based on the derivation and numerical examples presented above, the following conclusions can be drawn:(1)The existence of the other tunnel has influences on the vibration caused by the moving subway. And when the twin tunnels are both loaded, the dynamic responses cannot be calculated by simply summing up the results of a single tunnel model. Therefore, it is necessary to consider a twin-tunnel model when predicting the dynamic responses induced by subways for all load velocities. The dynamic responses of the single-phase ground are larger than those in the saturated ground which denotes that using a single phase model to predict the vibrations is more conservative.(2)The influence of the spacing between two tunnels on the dynamic response of the tunnel without the load applied is greater than that on the ground surface. With the increase of the spacing , the dynamic response of the adjacent tunnel decreases. Tunnel buried depth has a great influence on the dynamic response of adjacent tunnel and the ground surface. For points on the ground surface, the dynamic response decreases with the increase of the tunnel buried depth. Due to the presence of surface reflection region, the dynamic response of adjacent tunnel increases at first and then decreases as the burial depth increases. So, more attention should be paid to the spacing between tunnels and the tunnel buried depth in the construction of subway.(3)The subways moving at a low speed will generate a smaller dynamic response in the saturated half-space than moving at a higher speed. And when the load velocity approaches the Rayleigh wave velocity of the ground, which is around , the responses of the soil and liner will be amplified. Therefore, the subway moving speed is also an important parameter that affects the environmental vibration.
Appendix
The expressions of the elements in the mass matrix, the damping matrix, the stiffness matrix, and the force vector are presented as follows:where and , is the element interpolation function:
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work is supported by the National Key R&D Program of China (Grant No. 2016YFC0800200), the Projects of International Cooperation and Exchanges NSFC (Grant No. 51620105008), National Science Foundation for Young Scientists of China (Grant No. 51608482), and National Natural Science Foundation of China (Grant No. 51478424).