Abstract
Seismic evaluation of underground structures such as tunnels requires nonlinear dynamic analysis, due to the complex dynamic behavior of soil and the interaction of soil and structure. Simulation of the seismic response of the structure using nonlinear dynamic analysis is possible only with proper acceleration time history. Considering the vertical component of the earthquake (such as near-fault earthquakes) on the site is an important factor to achieve real structural responses. In the current study, soil-tunnel system has been modeled in ABAQUS software, considering Mohr–Coulomb nonlinear model for soil and concrete damage plasticity model for tunnel lining. In order to investigate the effect of seismic components correlation under different combinations of loads on the acceleration, axial force, and maximum shear force in tunnel lining, nonlinear dynamic analysis has been performed under four near-field earthquakes with different horizontal and vertical component ratios, considering 15 load combinations. The results show that increasing the vertical-horizontal component ratio has an insignificant effect on the maximum horizontal acceleration experienced by the tunnel lining. Also, the results of axial forces and shear forces indicate that increasing the ratio of vertical to horizontal components of the earthquake is the most effective factor on the axial force response.
1. Introduction
Due to the growing demand of underground structures, the high construction cost of them, their importance in the intermunicipal and interurban transport network, and the risk of injury caused by their damage, it is essential to study their stability against the earthquake. Although studies indicate that underground structures are more resistance than ground-based buildings against earthquakes, severe earthquakes can lead to catastrophic damage to underground structures as part of vital facilities in the countries [1]. Damage observed in the earthquakes of Lomaperita 1989, Kobe 1995, Duzce 1999, Chi Chi 1999, Niigata 2004, Wenchuan 2008, and Tohoku 2011 confirmed that underground structures are vulnerable to earthquakes [2]. Wang works in 1993 are considered to be one of the first steps in classifying past reports on the seismic performance of underground structures [3]. Few studies have been conducted to investigate the seismic behavior of tunnels constructed in rock or dense noncohesion soil [4]. Shallow tunnels constructed in softer soil are more vulnerable compared to deep tunnels built in rock [5, 6]. Underground structures have distinctive characteristics in comparison with on-ground structures. Being fully buried in the soil or rock and their significant length compared to the other dimensions of the tunnel are among the most distinctive features. As a result, the design of underground structures against seismic loads requires different considerations compared to the seismic design of ground structures. In underground structures, such as tunnels with a significant degree of containment due to earth conservation, contrary to on ground structures that appear as earthquake forces in the form of inertial forces, design and analysis is based on an attitude stating that its main criterion is the deformation of the earth and structure. These factors contribute to the better seismic performance of underground structures compare to ground structures [7].
Near the active faults, ground change is strongly influenced by the fault mechanism due to fault rupture with respect to the site (for instance, forward directivity and static permanent deformation at the fault location, which is known as filing effects). One of the most important distinctions of near-field earthquakes in comparison with conventional records is the long-range dynamic pulse with lasting ground deformation. The effects of the forward directivity are the main factor in the formation of the dynamic motion of the pulse with a large period. This pulse-like nature of the propagation is usually formed in the horizontal direction, and in the vertical direction of the fault. The effect of forwarding directivity often causes the component perpendicular to the fault to be larger than the parallel component of the fault, which is the case for the period of frequencies greater than 0.5 Hz. So far, extensive studies have been carried out to define the near-field mapping and its differentiation with a distant earthquake. In one category, the near-field mapping is distinguished based on engineering judgment from earthquakes far from the fault. It can be detected easily, especially if the earthquake speed mapping is available. Another criterion for detecting a near-field earthquake is the site distance to the seismic source. Usually, in near-field earthquakes, 15–30 km distances are defined as near-field. According to studies by Baker, a general definition for the near-field earthquake is presented. Based on this definition, the three points should be considered simultaneously to assign a near-field earthquake. These criteria include [8] the pulse index be greater than 0.85, pulses be formed at the early moments of the mapping speed, and the PGV earthquake record be larger than 30 cps. The vertical component of the earthquake is one of the most important parameters in near-field earthquakes. The effective parameters of the vertical component of the earthquake are as follows:(i)The earthquake magnitude: the vertical component increases with the rise in the magnitude of the earthquake response spectrum(ii)Distance from the fault: by increasing the distance from the fault, the value of the vertical component spectrum decreases(iii)Type of the fault: the vertical and horizontal spectrum values for reverse faults with superficial and high focal depth are larger than the values for the slip fault, and this corresponds to each other for periods larger than 1 second for different faults [9]
The vertical component of the earthquake was often ignored in the structural design; but slowly, with increasing acceleration of near-fault mappings and observing the large ratio of acceleration to the horizontal in near-fault earthquakes and observing empirically the destructive effects of vertical vibrations, this trend was changed and the vertical component effect was studied in the design of structures. One of the objectives of the tunnel cross-section selection is to distribute the uniform stresses in the soil mass around the tunnel. This issue has been widely discussed in technical resources. The drilling cross-sections that lead to the uniform distribution of tension around the tunnel are known as “harmonic holes.” More specifically, a drilling pattern in which the compressive stress across the roof and the walls are the same creates an optimal condition for the stress field. In a circular section in terms of the stress field, the distribution of compressive stress is often the same. Therefore, the greater the arc and curvature created in the tunnel cross-section, the tunnel’s seismic resistance increases [10], considering the nonlinear behavior of the soil and tunnel reduces the error in estimating seismic responses of the tunnel. Asgarian et al. researched performance evaluation of different types of steel moment resisting frames subjected to strong ground motion through incremental dynamic analysis and explained the effect of different types of ground motions on the results [11]. Ghasemzadeh and Abounouri researched compressional and shear wave intrinsic attenuation and velocity in partially saturated soils and studied the behavior of compressional (P1-wave) and shear wave velocities and intrinsic attenuations in sand saturated by two immiscible fluids (air and water) for a wide range of frequencies [12].
If the created deformation is the only effect of passing the seismic waves, from a structural point of view, flexibility is considered as a requirement for the design of the tunnel structure.
The use of nonlinear dynamic analysis in simulating the actual seismic response of a structure is possible only due to the time history of acceleration suitable for local soil conditions. The dynamic loading procedure is one of the important points to consider in dynamic analysis. This loading is usually possible by applying an earthquake to the bedrock level. In this paper, using the ABAQUS finite element software, Tehran’s soil and underground tunnel system are modeled using nonlinear time history analysis of the near-field fault records. The innovation of this research is taking into account different proportions of horizontal and vertical seismic components and the effect of correlation of seismic components of tunnel lining.
2. Seismic Analysis and Design of Underground Structures
When the earth deforms due to the passage of seismic waves, all underground structures also show a corresponding deformation. If the created deformation is the only effect of passing the seismic waves, from a structural point of view, flexibility is considered as a requirement for the design of tunnel structures. Therefore, in designing the correct and efficient structure of the tunnel, in addition to the strength capacity of the structural members, the structure and its flexibility should be considered. Based on Owen and Scholl studies, in underground structures, seismic vibrations are considered to hold three types of deformations, including axial compression and extension deformations, longitudinal bending, and ovaling/racking deformation [13]. When the earthquake wave propagated perpendicular to the tunnel axis, the stresses result in the formation of shear deformations of the type racking and ovaling in the tunnel cross-section. Design requirements in this type of deformation are usually applied to the width of the tunnel. The general behavior of the lining can be simulated as a buried structure subject to ground deformations under a two-dimensional plane strain condition. In general, the transmission waves cause bending in the structure, which creates alternate regions of compression and tension along the tunnel. The beam-like structure of the tunnel lining will then experience tension and compression on opposite sides. The analysis and design of underground structures are based on the structure deformations and the earth because the earth and the tunnel structure response is highly sensitive to earthquake-induced deformations. The methods of earthquake analysis of underground structures are experimental methods, physical modeling method, mathematical methods, quasi-static analysis, and numerical methods. Numerical methods, which are performed using lumped mass methods, finite difference, finite element, etc., are divided into two categories: pseudodynamic analysis and dynamic analysis. Another method for analyzing tunnels and underground structures is the dynamic time history analysis [1]. In the dynamic time history analysis, the entire system of structure and soil is exposed to dynamic stimulation in the soil-structure system. Also, the boundary conditions of the model should be chosen in such a way that the absorber of the seismic energy is induced and does not reflect the seismic wave again on the structure. Wang [3] and Penzien [14] proposed closed-form analytical solutions to calculate the thrust and bending moment in the tunnel lining and the effect of racking and ovaling. Hashash et al. [1, 15] and Patil et al. [16] explored the conventional methods for determining the seismic forces in the design of tunnels, taking into account different amounts of soil and tunnel interaction. In recent years, lack of knowledge in understanding the seismic behavior of the tunnel has led researchers such as Lanzano et al. [17] and Tsinidis [18, 19] to carry out research based on centrifuge model tests. Kalantarian and Dehghani studied the dynamic vibration relationship in the metro tunneling and proposed the propagation of Rayleigh waves as an input wave for issues in the time domain [20]. They concluded that, in horseshoe tunnel, distortion of cross-section as a combination of circular and rectangular section deformations occur. Also, the ground motion in the vertical direction is generally much less than its horizontal direction. Typically, vertical ground motion parameters are assumed to be 0.5 to 0.67 horizontal modes [21]. Abdel-Motaal and El-Nahhas studied the interaction between the tunnels and the surrounding granular soil and observed that the maximum strain in the tunnel lining was directly related to stiffness between the tunnel and surrounding soil (lining thickness and soil shear modulus) as well as the peak ground acceleration and the tunnel location (embedment depth). They conclude that a seismic analysis should be carried in regions subjected to peak ground acceleration greater than 0.15 g [13].
3. Numerical Modeling
The verification of the modeling was carried out in the ABAQUS [22] finite element software using Kouretzis et al. research, the analysis of circular tunnels due to seismic P-wave propagation, with emphasis on unreinforced concrete liners [23]. Kouretzis et al. modeled the tunnel and part of the surrounding area and analyzed the whole complex. In Kouretzis et al. research, the rock mass was 40 × 40 m2 and a circular section with a diameter of 8 m in its center was used. Rock and tunnel modeling is carried out in two-dimensional and in-plane strain conditions. The thickness of the tunnel lining is 15 cm with unreinforced concrete and is modeled with a wire element in the ABAQUS software. Elastic 2D beam elements (type = B21) were used with 64 elements in the perimeter to model tunnel lining. 2D plane strain continuum elements (type = CPE4) were used to model the surrounding rock and infinite elements (type = CINPE). Due to the half-infinity of the real model, the seismic energy applied in the model must go through the boundaries. In this case, it is necessary to encircle the surrounding using the energy absorbing boundaries to prevent the reflection of the waves in the area. Researchers have proposed different boundaries for dynamic analysis. Kouretzis et al. used Lysmer absorbing boundaries to model the energy absorbing boundaries up and down the model (20 meters from each side). Note that no loads from rock mass relaxation are assumed to be transferred to the final lining so that all the rock mass loads are borne by the temporary support shell. This assumption is appropriate for an unreinforced tunnel section constructed with the NATM method in the competent rock mass, where the final lining is installed after the primary lining has reached equilibrium [23]. The prediction of contact behavior between tunnel lining and surrounding is one of the key issues in the simulation of the tunnel response. Many linear elastic analytical solutions assume one of two zero friction modes (full-slip conditions) or a complete connection between the tunnel and the surrounding soil (no-slip conditions). Hence, it cannot correctly simulate the response of the contact surface to periodic loads. According to Wang, for most tunnels, contact conditions are a state between full slip and no slip because the full-slip condition reduces the axial force and the calculated forces under no-slip conditions may be several times greater than the full-slip rate [15]. Part of Kouretzis et al.’s research simulated no-slip conditions at the rock mass-structure interface, where the beam elements representing the liner are tied to the nodes of the surrounding rock mass. Elastic properties of the rock mass, tunnel lining, and seismic excitation are shown in Table 1. Tunnel lining geometry and finite element model of rock mass and tunnel in ABAQUS software are presented in Figure 1.
After performing the verification in no-slip condition and determining the maximum bending moment in tunnel lining and then comparing it with the proposed relationship provided by Kouretzis et al. [23] in Table 2, Table 3 shows that the maximum error generated between the responses is at an acceptable level of 10%.
Where
In another section, the study of Singh et al. was used for the verification of modeling and time history analysis [24]. They studied the seismic response in Delhi’s subway tunnels by finite element Plaxis software and Rayleigh’s damping based on horizontal and vertical components in Uttarkashi’s earthquake. The underground tunnels have been excavated through alluvium deposits, generally known as Delhi silt, and variation of elastic modulus of Delhi Silt with depth are shown in Figure 2 and Table 4 [24]. In-situ unit weight, γbulk, and saturated unit weight, γsat, were found to be 18 and 20 kN/m3, respectively. There are no underground water levels in the studied soil. Cohesion, c, of Delhi silt has been taken as zero (cohesionless silt). Friction angle, ϕ, and dilatational angle, ψ, were found to be 35° and 5°, respectively. The soil damping matrix is required for soil and structural analysis. According to the Rayleigh method and considering damping of soil as 5%, the coefficients α and β of the damping matrix of the studied soil are 0.161 and 0.013, respectively. The soil and tunnel modeling was carried out in two-dimensional and plane strain conditions, in which a rock mass of 140 × 60 m2 and a circular section with a diameter of 6.26 m at a depth of 20 m were used. The reinforced concrete tunnel lining thickness is 28 cm, its elastic modulus, Ec, is 3.16 × 107 kPa, Poisson’s ratio is 0.15, and the coefficients α and β of the damping matrix of the reinforced concrete lining of the tunnel with a 2% damping, according to the Rayleigh method, is 0.064 and 0.005, respectively.
Results are shown in Figures 3 and 4. As shown in Figures 3 and 4, the diagrams are close to each other.
In this paper, using the ABAQUS finite element software Teharn Metro Underground Tunnels and soil system using time history analysis under four near-field fault records of Tabas, Northridge, Lomaperieta, and Cape earthquakes, taking into account different ratios of vertical to horizontal components (αH + βV) modeled in 15 load combinations for each record. The underground tunnels which have been excavated through rock mass, generally known as Tehran soil and variation of elastic modulus of Tehran soil with depth, are shown in Figure 5 and Table 5. In-situ unit weight, γbulk, is 21.1 kN/m3. There are no underground water levels in the studied soil. Cohesion, c, of Tehran soil has been taken as 0.04 for GC and 0.07 for SC. Friction angle, ϕ, and dilatational angle, ψ, were found to be 33° and 3° for GC and 32° and 3° for SC, respectively. The soil damping matrix is required for soil and structural analysis. For the present case study, Rayleigh damping has been calculated. According to the Rayleigh method and considering damping of soil as 5%, the coefficients α and β of the damping matrix of the studied soil are 0.062 and 0.037, respectively. The soil and tunnel modeling was carried out in two-dimensional and plane strain conditions, in which a rock mass of 180 × 40 m2 and a circular section with a diameter of 8.85 m at a depth of 23.575 m were used. The reinforced concrete tunnel lining thickness is 35 cm, its elastic modulus, Ec, is 3.08 × 107 kPa, Poisson’s ratio is 0.2, and the coefficients α and β of the damping matrix of the reinforced concrete lining of the tunnel with 2% damping, according to the Rayleigh method, is 0.064 and 0.005, respectively, which is modeled with wire element using ABAQUS software. Elastic 2D beam elements (type = B21) were used in the perimeter to model tunnel lining. 2D plane strain continuum elements (type = CPE4R) were used to model the surrounding soils. Infinite elements (type = CINPE4), which are implemented in ABAQUS/standard to function as absorbing boundaries in wave propagation problems, were introduced at the left and right boundaries of the model, to avoid wave reflection phenomena and used in Lysmer absorbing boundaries to model the energy absorbing boundaries on the left and right sides of the model (70 meters from each side), and no-slip conditions at the soil mass-structure interface were simulated and the beam elements representing the liner are tied to the nodes of the surrounding soil mass.
To determine the nonlinear characteristics of tunnel lining in ABAQUS software, concrete damage plasticity behavioral model was used; for soil modeling, Mohr–Coulomb behavioral model was used. Table 6 shows the nonlinear tensile properties of concrete, and the specifications of concrete plastics are given in Table 7 [25].
The characteristics of the earthquake record are given in Table 8, and loading combinations in the three distinct groups with a constant coefficient of the horizontal component and variable coefficient of the vertical component are given in Tables 9–11.
Default parameters of tunnel lining and soil in the finite element model are given in Table 12.
The flowchart for research methodology is shown in Figure 6:
In the following, the effect of correlation of seismic components was investigated under different combinations of loads in determining maximum acceleration and yield stress, according to von Mises yield criterion, strain, axial force, shear force, and bending moment in tunnel lining.
4. Evaluation of Analytical Models
4.1. Effects of Load Combinations on the Maximum Acceleration Applied to the Tunnel Lining
In Figures 7–9, a comparison between the maximum vertical acceleration on different sections’ tunnel lining was performed under three near-field records of Tabas, Northridge, Lomaperieta, and Cape earthquakes in the 15 loading combinations for each record.
Examining the gradient values of the line showed that the maximum vertical acceleration of the tunnel lining is in Figures 7–9. Under Lomaperieta and Cape records, increasing the vertical to the horizontal component ratio, increases the slope of the line, and the maximum gradient of the line belongs to the Lomaperieta earthquake in the load combinations of 11 to 15. As the PVA (peak vertical acceleration) of the Lomaperieta earthquake record is 0. 896 and of the Cape earthquake record is 0.739. Vertical acceleration values at the 50% level (average acceleration level) are introduced on tunnel lining in Table 13. Table 13 shows that, in load combinations of 1 to 5, for β/α of 0.75, the difference between the vertical acceleration on tunnel lining in the two Lomaperieta and Cape records is equal to 0.03 g, the difference in vertical acceleration on the tunnel lining in both Cape and Northridge records is equal to 0.33 g, and the difference in vertical acceleration on the tunnel lining in both Northridge and Tabas records is equal to 0.11 g. In load combinations of 6 to 10, for β/α equal to 1.0, the difference in vertical acceleration on tunnel lining in two Lomaperieta and Cape records is equal to 0.03 g, the difference in vertical acceleration on tunnel lining in the two recordings of Cape and Northridge is equal to 0.31 g, and the difference in vertical acceleration on tunnel lining in the two recordings of Tabas and Northridge is equal to 0.11 g.
In load combinations of 11 to 15, for β/α equal to 0.6, the difference in vertical acceleration on tunnel lining in two Lomaperieta and Cape records is equal to 0.02 g, the difference in vertical acceleration on tunnel lining in the two recordings of Cape and Northridge is equal to 0.03 g, and the difference in vertical acceleration on tunnel lining in the two recordings of TABAS and Northridge is equal to 0.14 g.
In Figures 10–12, comparison between the maximum horizontal acceleration on tunnel lining was performed under four near-field records of Tabas, Northridge, Lomaperieta, and Cape in 15 load combinations for each record. Investigating the gradient values of the horizontal gradient line in Figures 10–12 shows that the slope of the line has been slowly reduced or increased by increasing the vertical to the horizontal component ratio. This indicates the low impact of the vertical components in the maximum horizontal acceleration on tunnel lining, and the highest gradient of the line belongs to the Lomaperieta and Cape earthquakes. The reason is that the PHA (Peak Horizontal Acceleration) of the Lomaperieta earthquake record is 0.607 and of the Cape earthquake record is 1.493.
4.2. The Effect of Load Combinations on the Forces in Tunnel Lining
The effects of change in the ratio of seismic components on the maximum axial force of tunnel lining are shown in Figures 13–15.
The review of Figures 13–15 shows that when the Cape earthquake is applied to the soil and tunnel set, the highest axial force values belong to the load combinations of 12 to 14, and the maximum F line slope is −262.5 and belongs to the load combinations of 11 to 15. The review of Figures 13–15 shows that when the Lomaperieta earthquake is applied to the soil and tunnel set, the highest axial force values belong to the load combinations of 11 to 15, and the maximum F line slope is 860 and belongs to the load combinations of 11 to 15. The review of Figures 13–15 shows that when the Northridge earthquake is applied to the soil and tunnel set, the highest axial force values belong to the load combinations of 11 to 15 and the reason for this is a horizontal component earthquake coefficient of 1.25 that indicates the importance of the horizontal component of the earthquake in determining the value of axial force tunnel lining. The maximum F slope is −84 and belongs to the loading combination of 11 to 15. By comparing the effects of Tabas, Northridge, Lomaperieta, and Cape earthquakes on the soil and tunnel set shows that the greatest axial force created in the tunnel lining belongs to the Cape earthquake, due to the larger PVA and PHA of the cape earthquake than the other earthquakes.
Different seismic waves have different Fourier spectra and different predominant frequencies. This effect on the dynamic response of tunnel structure cannot be ignored. When “β/α” is the same, the influence of seismic waves on the dynamic response of the tunnel is different. The review of Figures 13–15 shows that the maximum F line slope belongs to the load combinations of 6 to 10 (α = 0.75 and β = Var).
The effects of change in the ratio of seismic components on the maximum shear force of tunnel lining (V) are shown in Figures 16–18.
The review of Figures 16–18 shows that when the Tabas earthquake is applied to the soil and tunnel set, the highest shear force values belong to the load combinations of 11 to 15, and the maximum line slope is 5.88 and belongs to the load combination of 6 to 10. The review of Figures 16–18 shows that when the Northridge earthquake is applied to the soil and tunnel set, the highest shear force values belong to the load combinations of 11 to 15. The line slope in the load combinations of 1 to 15 is decreasing. The highest shear force values belong to the load combinations of 11 to 15.
The review of Figures 16–18 shows that when the Cape earthquake is applied to the soil and tunnel set, the highest shear force values belong to the load combinations of 11 to 15. The line slope of the shear force in the loading combination 11 to 15 is decreasing. the maximum line slope is 65 and belongs to the load combinations of 11 to 15. As seen in Figures 16–18, increasing the ratio (β/α) leads to a decrease or increase the , depending on the earthquake type and three load groups.
This is due to the difference between the characteristics of the three various regions in the time history of the strong ground motion in the considered accelerators because all the effective seismic ground motion signals consist of three regions with completely different characteristics in which the first region has shown to be mild with increasing slope, the second region generally has much constant energy and power, and the third region has decreasing circumstance. Therefore, with respect to the variety of time classification in separating the three regions among the selected records in the dynamic analysis of tunnels, the alteration regime in the force components has different conditions depending on the coefficient of the vertical to the horizontal component [26].
5. Conclusion
The effect of the earthquake on the soil and tunnel set depends on various parameters such as the maximum horizontal and vertical acceleration of the earthquake, the magnitude and duration of the earthquake, and the relative stiffness between the tunnel and the ground and tunnel lining. Since the vertical component of the earthquake is one of the most important parameters in the near-field fault earthquake, in this study, using nonlinear time history analysis under four near-field recordings considering different horizontal and vertical component ratios in 15 loading combinations, the seismic behavior of the tunnel with reinforced concrete lining surrounded through rock mass (Tehran soil) was performed:(i)The PVA value is the most important factor in increasing the maximum vertical acceleration in tunnel lining in all load combinations, and since the Lomaperieta earthquake has the highest PVA, it has the highest vertical acceleration in all load combinations(ii)The PHA is the most important factor in increasing the maximum horizontal acceleration in tunnel lining in all load combinations, and PVA has little effect, since the Cape earthquake has the highest PHA, thus having the highest horizontal acceleration in all load combinations(iii)With increasing vertical component in four separate load groups, the maximum axial force value has increased and the highest increase in the axial force is due to the application of the Cape earthquake to the soil and tunnel set and in the load combinations of 11 to 15(iv)The increase of the vertical component in three separate load groups will lead to a decrease or increase in the value of shear force, and the largest increase in shear force is due to the application of the Cape earthquake to the soil and tunnel set and in the load combinations of 11 to 15
Data Availability
The Tehran metro underground tunnel and soil system data used to support the findings of this study have been deposited in the TTBP Consulting Engineers Company repository (TTBP Consulting Engineers Company (2015). Design of T-Junction Subway Line Connector Structure; Case Study: Line 3 and 4 in Tehran Metro1. Quarterly Journal of Boland Payeh Company, Vol. 3, No. 10). Link: bolandpayeh.com › Journals › BolandpayehJ10.
Conflicts of Interest
The authors declare that they have no conflicts of interest.