On subsurface box-shaped lined tunnel under incident SH-wave propagation

Abstract

In this paper, a half-plane time-domain boundary element method is applied to obtain the seismic ground response, including a subsurface box-shaped lined tunnel deployed in a linear homogenous elastic medium exposed to obliquely incident SH-waves. Only the boundary around the tunnel is required to be discretized. To prepare an appropriate model by quadratic elements, a double-node procedure is used to receive dual boundary fields at corners as well as change the direction of the normal vector. After encoding the method in a previously confirmed computer program, a numerical study is carried out to sensitize some effective parameters, including frequency content and incident wave angle for obtaining a surface response. The depth and impedance ratio of the lining are assumed to be unvaried. The responses are illustrated in the time and frequency domains as two/three-dimensional graphs. The results showed that subsurface openings with sharp corners distorted the propagation path of the anti-plane waves to achieve the critical states on the ground surface. The present approach can be proposed to civil engineers for preparing simple underground box-shaped models with angular boundaries.

Graphical abstract

Keywords

box-lined tunnel / half-plane BEM / surface response / SH-wave / time-domain

1 Introduction

The local site effect is considered as one of the most decisive factors in earthquake occurrence, and its importance and complexity lead researchers to study more comprehensively. The development of the transport industry and the pattern of major cities have increased the importance of underground tunnels such as transportation tunnels. Thus, the behavior knowledge of the ground surface in the existence of subsurface openings is one of the preliminary studies for constructing any surface structures. In this regard, simple and accurate modeling of such subsurface effects is considered an urgent need. Various methods have been referred for seismic analysis of topographic features and determination of ground surface response. In the literature, some empirical approaches proposed by Sgarlato et al. [1] and Kazemeini et al. [2] are observed for investigating the site effects in the presence of underground tunnels. In general, seismic analysis methods of topographic features have been separated into analytical, experimental, and numerical methods [3]. Although in the analytical methods, the responses are more accurate, the models are limited to a series of simple geometries. First, Asano [4] and then Datta [5] studied the effects of topography under SH-wave using analytical methods. With the development of the wave function expansion and weighted residual methods, the analytical methods were used in the analysis of underground tunnels exposed to SH-waves by the majority of researchers, including Gamer [6], Lee and Trifunac [7], Chen [8], Datta and Shah [9], Lee and Manoogian [10], Shi et al. [11], Manoogian [12], Jiang et al. [13], Liang et al. [14], Yi et al. [15], Amornwongpaibun et al. [16], Kara [17] and Gao et al. [18]. Another analytical procedure used by researchers (e.g., Wang and Liu [19] and Shi et al. [20]) was a multi-polar coordinate method. In the meantime, analytical studies of Lee et al. [21] and Smerzini et al. [22] were performed by Bessel function expansion and Liu and Lin [23] and Tsaur and Chang [24] used the region-matching technique. Recently, He et al. [25] analytically modeled the vibration prediction of tunnels buried in a full space. Moreover, the volumetric numerical methods (i.e., domain approaches) have been developed in the static/seismic analysis of ground surface, including subsurface tunnels. The technical literature shows that Lee and Rowe [26], Molinero et al. [27], Karakus et al. [28], Yamamoto et al. [29], Zhang et al. [30], Das et al. [31], and Tsinidis [32] were among the researchers who studied the tunnel problems by FEM. The use of FDM can name the studies of Yiouta-Mitra et al. [33], Chakeri et al. [34], Besharat et al. [35], Do et al. [36], Jamshasb and Shahrbabaki [37], Narayan et al. [38], and Ziaei and Ahangari [39] as well. Although the domain approaches have abilities including appropriate accuracy, perceptible mathematics, and development on various issues, they have time-consuming calculations as well as complicated sketches for modeling the problems with unlimited boundaries. Therefore, these issues paved the way for the emergence of boundary approaches.

In the boundary element method (BEM), only boundaries around the considered features need to be discretized. It is an appropriate method for dynamic analysis of geotechnical problems because the wave’s radiation conditions are automatically satisfied in the formulation [40]. BEM approaches have always been used for modeling the subsurface cavities and underground tunnels subjected to seismic waves as well as static loads. Luco and de Barros [41], Yu and Dravinski [42], Parvanova et al. [43], and Liu and Liu [44] were among the researchers who used full-plane frequency-domain BEM to examine the effect of subsurface cavities on the response of the ground surface. Benites et al. [45] determined the response of the ground surface, including multiple cavities with a frequency-domain half-plane BEM. In the use of elastostatic full-plane/half-plane BEM, Panji et al. [46,47], and Panji and Ansari [48] presented the response of deep and shallow tunnels subjected to effective loads. Because of advantages such as nonlinear analysis of different problems in combination with other numerical methods, obtaining real solutions, and analysis of time-dependent geometries, BEM was developed in the time-domain. The full-plane time-domain BEM was used by Takemiya and Fujiwara [49] and Kamalian et al. [50] in analyzing the effects of surface topographies. Also, the half-plane time-domain BEM was proposed by some researchers, including Rice and Sadd [51], Belytschko and Chang [52], Hirai [53], and used by Panji et al. [54,55,56] for analyzing seismic ground motion affected by subsurface topographic features and underground tunnels. Recently, Panji and Ansari [57] and Panji and Mojtabazadeh-Hasanlouei [58,59] studied the SH-wave scattering by the lined and unlined tunnels located in an elastic half-plane, respectively. The scattering of plane P and SV-waves by lined tunnels embedded in an elastic half-space was discovered by Huang et al. [60]. In the most recent studies, Panji et al. [61] and Mojtabazadeh-Hasanlouei et al. [62] studied the SH-wave dispersion by single and multiple inclusions, respectively. Also, the SH-wave propagation by periodically distributed semi-sine-shaped canyons was investigated by Panji and Mojtabazadeh-Hasanlouei [63]. Among the recent works discussing the wave propagation problems with BEMcan name the studies of Chen et al. [64,65] and Shaaban et al. [66,67] as well.

The literature review showed that the lined tunnels could be used for mining tunnels and underground transportation systems like Metros, etc. [68]. As far as the authors know, a general seismic investigation was not done on the box-shaped (square-shaped) lined tunnels placed in an elastic half-plane. Therefore, in this paper, the seismic ground motion is obtained for the mentioned tunnel exposed to oblique incident SH-wave propagation using half-plane time-domain BEM. The primary purpose of this paper is focused on obtaining an accurate response of angular-shaped subsurface structures and their simple modeling by time-domain half-plane BEM.

2 Problem statement

Figure 1 shows a box-shaped lined tunnel placed in a linear homogeneous elastic isotropic half-plane under incident out-of-plane SH-waves of the Ricker type. As can be observed, only discretizing the tunnel boundary was sufficient for modeling the tunnel with the assistance of half-plane time-domain BEM. In this figure, \Omega is the domain and the boundaries of the body are defined by \mathrm{\Gamma }. The dimensions of the square-shaped tunnel are shown by 2b, where b is half the length of each side. The wave velocity of the half-plane is presented by c. \theta is the angle of the incident waves and {\mathit{\lambda }}_{\mathrm{p}} is the wavelength. The input motion function is defined as the Ricker wavelet type presented in Eq. (1):

(1)f\left(t\right)=\left[\mathrm{1}-\mathrm{2}{\left(\mathrm{\pi }{f}_{\mathrm{p}}\left(t-t_{\mathrm{0}}\right)\right)}^{\mathrm{2}}\right]{\mathrm{e}}^{-{\left(\mathrm{\pi }{f}_{\mathrm{p}}\left(t-t_{\mathrm{0}}\right)\right)}^{\mathrm{2}}},

where the dominant frequency of input motion is presented by f_{\mathrm{p}} and t_{\mathrm{0}} is the time-transition factor. Figure 2 shows the Ricker wavelet in time and frequency-domain, respectively. Based on a half-plane model, the free-field displacement (u^{\mathrm{ff}}) can be obtained as follows [54,56]:

u^{\mathrm{ff}}\left(x,y,t\right)={\alpha }_{\max }

(2)·\left(\begin{array}{l}\left[\mathrm{1}-\mathrm{2}{\left({\displaystyle \frac{\mathrm{\pi }{f}_{\mathrm{p}}}{c_{\mathrm{1}}}}{\alpha }^{\mathrm{inc}\mathrm{.}}\right)}^{\mathrm{2}}\right]e^{-{\left(\frac{\mathrm{\pi }{f}_{\mathrm{p}}}{c_{\mathrm{1}}}{\alpha }^{\mathrm{inc}\mathrm{.}}\right)}^{\mathrm{2}}}H\left(t-{\displaystyle \frac{\left|y\right|}{c_{\mathrm{1}}}}\right)+\hfill \\ \left[\mathrm{1}-\mathrm{2}{\left({\displaystyle \frac{\mathrm{\pi }{f}_{\mathrm{p}}}{c_{\mathrm{1}}}}{\alpha }^{\mathrm{ref}\mathrm{.}}\right)}^{\mathrm{2}}\right]e^{-{\left(\frac{\mathrm{\pi }{f}_{\mathrm{p}}}{c_{\mathrm{1}}}{\alpha }^{\mathrm{ref}\mathrm{.}}\right)}^{\mathrm{2}}}H\left(t-{\displaystyle \frac{\left|y\right|}{c_{\mathrm{1}}}}\right)\hfill \end{array}\right),

where {\alpha }_{\max } is the maximum time-history amplitude, u^{\mathrm{ff}}\left(x,y,t\right) is the free-field displacement at the point (x,\mathrm{\hspace{0.17em}}z) and time t, H is the Heaviside function, {\alpha }^{\mathrm{inc}\mathrm{.}} and {\alpha }^{\mathrm{ref}\mathrm{.}} are the phase of incident and reflected waves as follows:

(3){\alpha }^{\mathrm{inc}\mathrm{.}}=c\left(t-t_{\mathrm{0}}\right)+\left|y\right|,

(4){\alpha }^{\mathrm{ref}\mathrm{.}}=c\left(t-t_{\mathrm{0}}\right)+\left|y\right|.

3 Half-plane time-domain BEM

Two-dimensional scalar wave equation and governing boundary conditions on the ground surface were respectively introduced for a linear homogeneous elastic medium as follows [69,70]:

(5){\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_j\left(x,y,t\right)}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_j\left(x,y,t\right)}{\partial y^{\mathrm{2}}}}+b_j\left(x,y,t\right)={\displaystyle \frac{\mathrm{1}}{c_j^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_j\left(x,y,t\right)}{\partial t^{\mathrm{2}}}},j=\mathrm{1},\mathrm{2}

(6){{\displaystyle \frac{\partial u_{\mathrm{1}}\left(x,y,t\right)}{\partial n}}|}_{y=\mathrm{0}}=\mathrm{0.}

In the above equation, c_j is the shear wave velocity, u_j\left(x,y,t\right) and b_j\left(x,y,t\right) are the antiplane displacement and body forces at coordinate \left(x,y\right) for time t and n is the surface normal vector. By considering boundary conditions in Eq. (6), the time-domain half-plane Green’s functions were achieved from a singular solution of Eq. (5) [54].

3.1 Boundary integral equation

By establishing the weighted residual integral to Eq. (5) and subsequently disregarding body terms and initial conditions, and considering principles of waves diffraction in a half-plane, the modified direct time-domain boundary integral equation (BIE) can be determined as follows [71,72]:

(7)c\left(\xi \right)u\left(\xi ,t\right)=\underset{{\mathrm{\Gamma }}_{\mathrm{1}}\cup {\mathrm{\Gamma }}_{\mathrm{2}}}{\overset{}{{\displaystyle \int }}}\{{{\displaystyle \int }}_{\mathrm{0}}^t[u^{*}\left(x,t\mathrm{;}\xi ,\tau \right)·q_j\left(x,t\right)-q^{*}\left(x,t\mathrm{;}\xi ,\tau \right)·u_j\left(x,t\right)]\mathrm{d}\tau \}\mathrm{d}\mathrm{\Gamma }\left(x\right)+u_{\mathrm{1}}^{\mathrm{ff}}\left(\xi ,t\right),j=\mathrm{1},\mathrm{2},

where u^{*}\left(x,t\mathrm{;}\xi ,\tau \right) is the transient half-plane displacement Green’s function in the position x and current time t due to an anti-plane single pulse at position \xi and time \tau. Also, q^{*}\left(x,t\mathrm{;}\xi ,\tau \right) is the transient half-plane traction Green’s function obtained from derivative of the displacement Green’s function, u_j and q_j are the boundary displacement and traction values, respectively, c\left(\xi \right) is the corner effect, and u_{\mathrm{1}}^{\mathrm{ff}} is the free field motion of the surface for a uniform and homogeneous half-plane, which also satisfies the stress-free boundary conditions. To solve Eq. (7), it is needed to discretize the time-axis and geometric boundary of the body. By dividing the time-axis into N equal parts and considering linear variations in each time interval, the temporal integration on Green’s functions can be analytically performed. By overlapping adjacent time-nodes and deleting singular terms of wavefront, the above BIE can be written as:

c\left(\xi \right)u^N\left(\xi \right)={{\displaystyle \sum }}_{n=\mathrm{1}}^N

\underset{{\mathrm{\Gamma }}_{\mathrm{1}}\cup {\mathrm{\Gamma }}_{\mathrm{2}}}{\overset{}{{\displaystyle \int }}}\left(\begin{array}{l}\left[U_{\mathrm{1}}^{N-n+\mathrm{1}}\left(x,\xi \right)+U_{\mathrm{2}}^{N-n}\left(x,\xi \right)\right]q_j^n\left(x\right)-\hfill \\ \left[Q_{\mathrm{1}}^{N-n+\mathrm{1}}\left(x,\xi \right)+Q_{\mathrm{2}}^{N-n}\left(x,\xi \right)\right]u_j^n\left(x\right)\hfill \end{array}\right)\mathrm{d}\mathrm{\Gamma }\left(x\right)

(8)+u_{\mathrm{1}}^{\mathrm{ff}.N}\left(\xi \right),j=\mathrm{1},\mathrm{2},

where U_{\mathrm{1}}^{N-n+\mathrm{1}}+U_{\mathrm{2}}^{N-n} and Q_{\mathrm{1}}^{N-n+\mathrm{1}}+Q_{\mathrm{2}}^{N-n} are the abbreviated closed-form kernels for half-plane time-domain displacement and traction of the scalar wave equation, respectively, which were proposed by Panji et al. [54]. Also, u^{\mathrm{ff}\mathrm{.}N} is the free-field motion of the surface in the time-step N, and u_j^n and q_j^n are the displacement and traction fields for the boundary nodes, respectively.

4 Numerical implementation

Before solving Eq. (8), the discretization of the geometric boundary of the body, as well as the time-axis, should be carried out for obtaining the field variables. Up to this point in solving the problem, the computational process was accurate. In this step, an analytical process and a numerical procedure were applied on Eq. (8) to obtain its matrix form for temporal and spatial integration, where the details of these processes can be found in Panji and Ansari [57,73].

4.1 Double-node procedure at corners

As depicted in Fig. 3, due to the sharp corners in the geometry of the box-lined tunnel, the double-node procedure has been used in the model preparation. Based on this process which is one of the capabilities of BEM, two independent nodes are considered at the sharp corners of the tunnel to form the independent columns of the corresponding matrix for each node. Using two independent nodes at sharp corners enables easy evaluation of the jump term (geometry coefficient) as all collocation points are placed on smooth surfaces. This process increases the accuracy of modeling and provides exact responses, which is previously used by Ji et al. [74] for TD-BEM and static study of Panji et al. [47].

5 Applications

The proposed method is implemented and developed in the general algorithm known as DASBEM [54]. This program is prepared for the seismic analysis of 2-D underground box-shaped lined tunnels embedded in an elastic half-plane. Also, the implementation of the numerical procedure has been done in the MATLAB [75] programming software. Some verification examples around the different shapes of lined tunnels were solved, which can be found in the studies of Panji and Ansari [57,73], and in the present work, this part was ignored. The time-domain responses of the ground surface were presented by some synthetic seismograms beside two and three-dimensional amplification patterns in the frequency domain. In this regard, incident wave angle and dimensionless frequency of the responses were studied as the intended parameters. Based on the following sections, it is necessary to introduce the utilized parameters. DR is the depth ratio of the tunnel and is fixed on 2.5 below the surface. Also, IR is the interspace ratio between the surface stations and equal to 2. The impedance ratio (I) is the rigidity ratio of lining to the adjacent space and defined as (I={\rho }_{\mathrm{l}}{c}_{\mathrm{l}}/{\rho }_{\mathrm{m}}{c}_{\mathrm{m}}) in which, {\rho }_{\mathrm{l}} and c_{\mathrm{l}} are the mass density and shear-wave velocity of lining; {\rho }_{\mathrm{m}} and c_{\mathrm{m}} are the mentioned parameters for the domain, respectively. It should be noted that the impedance ratio is considered equal to 10 and states that the concrete lining is 10 times harder than the surrounding medium. The other parameters of lining are applied according to the study of Panji and Ansari [73]. The amplification ratio is defined as the ratio of the surface response amplitude to free-field motion. The ratio of the Fourier amplitude of the total ground surface motion obtained by BEM to the Fourier amplitude of the incident motion is known as the normalized displacement amplitude (NDA). The dimensionless frequency (\eta) is determined by (\eta =\omega b/\mathrm{\pi }c), where \omega is the angular frequency of the wave, b is the width of the tunnel, and c is the shear-wave velocity of the domain. The dimensionless period (P) is defined as the inverse of dimensionless frequency (P=\mathrm{1}/\eta) as well.

In section 5.1, the 3-D and 2-D responses of the time-domain are demonstrating the diffraction of the incident waves on the ground surface in the presence of a box-shaped lined tunnel, respectively. Then, the 3-D and 2-D frequency-domain responses are illustrated the overall pattern of amplification and displacements of the surface in section 5.2. Moreover, the 2-D diagram of amplification versus dimensionless period is presented for different incident angles in the last part. One of the most effective factors to determine the acceptable responses is selecting the appropriate values for the nodal interval. To measure the convergence of responses, the sensitivity of the nodal interval was studied by \beta factor. \beta is defined as \beta =c\mathrm{\Delta }t/\mathrm{\Delta }l where, \mathrm{\Delta }t is the increment of time, and \mathrm{\Delta }l is the nodes interval. To calculate \beta factor, the values of 1200 m/s and 0.007 s are considered for c and \mathrm{\Delta }t for 1000 time-steps. By considering four different \mathrm{\Delta }l of 15, 25, 35, and 45 m, the values of 0.560, 0.336, 0.240, and 0.187 were obtained for \beta, respectively. Figure 4 shows the time-history response of the ground surface for different \beta due to vertically propagating incident SH-waves. It is observed that the responses are always stable by changing the distance of discretization on the inner and outer surfaces of the lined tunnel. The minimum nodal interval of 15 m is considered for all of the models. Although the analysis was carried out at other distances, the mentioned value is selected for the node-interval due to obtaining the convergent time-domain results as well as smooth responses at high frequencies. In the following, some of the responses are presented in the time and frequency domain.

5.1 Time-domain responses

Figure 5 shows the general diffraction of the waves on the ground surface in the time-domain exposed to obliquely propagating incident SH-wave. To obtain the synthetic seismograms, the impedance ratio for concrete lining is considered equal to 10 [73], the angles of 0°, 30°, 60°, and 90° are applied for the incident front, and the ground surface range is shown between −4b and 4b. The responses are displayed in 12 s. as well. For a better and accurate interpretation of the results, the different phases of incident waves are marked by some specific stations. In this regard, D shows the direct waves, R_L is the resonance of liner which is quickly propagated in whole parts of the liner, and its effect reaches the ground surface sooner than the wave-front. R illustrates the reflected phase due to the wave’s collision to the outer boundary of the liner. Also, C displays the crawler waves when they crept on the boundary of the tunnel. T corresponds to the trapped phase, where the waves experience intermittent reflections. Because of the smooth surfaces and hard material of the liner compared to the surrounding soils, it acts as a mirror and traps the transient waves between the boundary of the tunnel and surface. This phenomenon increases the fluctuations of the responses and extends the duration of convergence time. When the vertical incident front is applied to the model (Fig. 5(a)), in addition to the symmetry of response, the stations of R and R_L, and a little further, C and T are located in the same position, respectively. In this figure, the resonance effect of the liner is clearly observable where the oscillations are started before the collision of transient waves to the ground surface (location of station D). Continuation of vibrations after the second three is due to the trapped waves as well. When the incident angle is 30° (Fig. 5(b)), the effect of R_L is weak compared to the vertical case, and the path of crawler waves appears behind the tunnel. By increasing the inclination of the wavefront to 60° (Fig. 5(c)), a significant volume of seismic waves is reflected in the collision with the tunnel, and the stations of R_L, R, and T are concentrated. Also, the volume of crawler waves is weaker which is because of the specific shape of the tunnel. In horizontal radiation of waves (Fig. 5(d)), the effect of creeper waves on the surface is removed, and the most volume of waves is reflected directly when they reached the smooth surface of the tunnel. In the following, Fig. 6 shows the 2-D time-domain responses and the amplitudes of waves for three stations on the surface. The interspace range of stations one and three from the centerline is equal to 2. As can be seen in Fig. 6(a), the responses of the first and last stations are coincident with the amplitude of 2.2 for the vertical incident. By inclination of wavefront to 30° (Fig. 6(b)), the first station recorded the most amplitude with the value of 2, but the other stations are displayed a significant reduction of amplitude. This case clearly shows the mirror effect of the liner’s smooth surface in the reflection of seismic waves. By increasing the wavefront inclination to 60° (Fig. 6(c)), the amplitudes are the same as the previous case (Fig. 6(b)) with the maximum value of 2, but unlike the 30° case, the second station on the top of the tunnel also showed significant growth of amplitude. For the horizontal wavefront, the response of each station is high when compared to the previous cases, but the third station recorded the maximum amplitude with the value of 2.27 among the other results.

5.2 Frequency-domain responses

Although the time-domain results seem desirable to observe the dispersion of the wave, the only possible way to examining the amplification ratio and displaying its behavioral patterns is utilizing the frequency-domain responses. Thus, Fig. 7 shows the three-dimensional surface amplification pattern. The 3-D frequency-domain responses are presented for dimensionless frequencies of 0.25 to 8. As can be seen, when the incident front is applied vertically (Fig. 7(a)), the response pattern is symmetric, and the shadow zone [76] has appeared above the location of the tunnel. This phenomenon shows that the existence of the tunnel below the surface blocked the passage of seismic waves and created a safe area above the tunnel, but the highest seismic hazard has accumulated in the peripheral zones around the tunnel. In this case (Fig. 7(a)), the maximum amplification of 1.65 is recorded on the location of tunnel edges. By inclination of the wavefront to 30° (Fig. 7(b)), the shadow-zone is moved behind the tunnel, the turmoil of fluctuations is reduced on the opposite side of the wavefront, and the maximum amplification is increased to the value of 1.9 on the radiation side. When the incident front is applied at 60° (Fig. 7(c)), the effect of shadow-zone is stronger. Therefore, the accumulation of seismic waves has caused a higher amplification with a value of more than 2 on the propagation side. In 90° of wave propagation (Fig. 7(d)), the greatest amplifications emerge with the value of 3.4. Thus, in the horizontal wavefront, the tunnel has the opposite function and has led to a significant increase in surface amplification and the creation of a high-risk area on the location of the tunnel. To get a better view, Fig. 8 shows the 2-D ground surface displacement for dimensionless frequencies of 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.0, respectively. As shown in this figure, with increasing the dimensionless frequency, the number of fluctuations is increased. On the other hand, comparing the results shows that the maximum surface displacement of 6.2 occurs right above the location of the tunnel when the horizontal wave-front is applied to the model (Fig. 8(g)). But this value does not exceed 3.5 for other incident angles. This result is also clearly demonstrated the high seismic hazard on the location of the tunnel when exposed to horizontal wave propagation. In the last part, to determine the sensitivity of the amplification ratio in different incident angles, three stations are considered on the ground surface, and the response pattern is presented versus the dimensionless period. As illustrated in Fig. 9, according to the interest of engineering knowledge for studying the responses in the periodic range of 0.25 to 10 [77], which corresponds to the wavelengths of 0.25 to 10 times of the tunnel width, this range of dimensionless period is considered in the present study as well. According to the previous interpretations, the highest amplification of 2.1 is recorded in the second station when the wave-front was horizontal (Fig. 9(d)). The highest results occur in dimensionless periods of less than 1. Also, the responses are converged to the amplification of the smooth ground surface (unit value) at wavelengths larger than three times the tunnel width.

6 Conclusions

A uniform elastic homogenous half-plane medium with a box-shaped lined tunnel was analyzed to obtain the seismic response of the ground surface exposed to obliquely propagating incident SH-waves. A half-plane time-domain BEM previously proposed by Panji et al. [54] was developed to establish the underground-lined tunnel model. By focusing the meshes on the opening boundaries and avoiding the ground surface discretization, not only the calculation time was reduced, but also the volume of input data was minimized compared to traditional BEM approaches. Moreover, in the use of double-node procedure at corners, the box-shaped openings were modeled more realistic. Next, a numerical study was carried out using a lined box tunnel embedded in the assumed depth, which its results were summarized as follows.

1) The separated phases of transient waves were clearly showed in the 3-D synthetic seismograms for different wave-front angles.

2) The 2-D time-domain responses demonstrated that the maximum amplitude of 2.27 was recorded behind the location of the tunnel when the wavefront was a horizontal incident.

3) The overall pattern of the frequency-domain responses showed that, when the tunnel was exposed to the propagation of vertical SH-waves, the isolation effect of the tunnel was quite recognized on the reduction of the seismic ground motion.

4) When the angle of the wave’s front was equal to 60°, because of the shadow zone effect, the fluctuations were decreased behind the tunnel. This issue was weak in the horizontal incident waves.

5) The 3-D frequency-domain results showed that the maximum amplification ratio was recorded for horizontal wave-front with the value of 3.4.

6) The obtained results versus the dimensionless period showed that the highest amplification of 2.1 was achieved above the location of the tunnel in the horizontally radiating wave-front. Moreover, the responses were converged to the amplification of the smooth ground surface (unit value) at wavelengths larger than three times the tunnel width.

References

[1] Sgarlato G, Lombardo G, Rigano R. Evaluation of seismic site response nearby underground cavities using earthquake and ambient noise recordings: A case study in Catania area, Italy. Engineering Geology, 2011, 122(3–4): 281–291

[2] Kazemeini M J, Haghshenas E, Kamalian M. Experimental evaluation of seismic site response nearby underground cavities (study of subway tunnel in city of Karaj, Iran). Iranian Journal of Science and Technology. Transaction of Civil Engineering, 2015, 39: 319–332

[3] Sánchez-Sesma F J, Palencia V J, Luzón F. Estimation of local site effects during earthquakes: An overview. ISET Journal of Earthquake Technology, 2002, 39(3): 167–193

[4] Asano S. Reflection and refraction of elastic waves at a corrugated boundary surface. Part I, The case of incidence of SH-wave. Bulletin Earthquake Research Institute, 1960, 38: 177–197

[5] Datta S K. Diffraction of SH-waves by an elliptic elastic cylinder. International Journal of Solids and Structures, 1974, 10(1): 123–133

[6] Gamer U. Dynamic stress concentration in an elastic half space with a semi-circular cavity excited by SH-waves. International Journal of Solids and Structures, 1977, 13(7): 675–681

[7] Lee V W, Trifunac M D. Response of tunnels to incident SH-waves. Journal of the Engineering Mechanics Division, 1979, 105(4): 643–659

[8] Chen Y L. The analysis of elastic liner in a cylindrical tunnel subjected to SH-waves. Journal of the Chinese Institute of Chemical Engineers, 1980, 3(1): 21–29

[9] Datta S K, Shah A H. Scattering of SH-waves by embedded cavities. Wave Motion, 1982, 4(3): 265–283

[10] Lee V W, Manoogian M E. Surface motion above an arbitrary shape underground cavity for incident SH-waves. Europ Earthquake Engineering, 1995, 8(1): 3–11

[11] Shi S, Han F, Wang Z, Liu D. The interaction of plane SH-waves and non-circular cavity surfaced with lining in anisotropic media. Applied Mathematics and Mechanics, 1996, 17(9): 855–867

[12] Manoogian M E. Scattering and diffraction of SH-waves above an arbitrarily shaped tunnel. ISET Journal of Earthquake Technology, 2000, 37(1–3): 11–26

[13] Jiang L F, Zhou X L, Wang J H. Scattering of a plane wave by a lined cylindrical cavity in a poroelastic half-space. Computers and Geotechnics, 2009, 36(5): 773–786

[14] Liang J, Luo H, Lee V W. Diffraction of plane SH-waves by a semi-circular cavity in half-space. Earth Science, 2010, 23(1): 5–12

[15] Yi C, Zhang P, Johansson D, Nyberg U. Dynamic analysis for a circular lined tunnel with an imperfectly bonded interface impacted by plane SH-waves.‏ In: Proceedings of the World Tunnel Congress-Tunnels for a better Life. Foz do Iguaçu: World Tunnel Congress, 2014

[16] Amornwongpaibun A, Luo H, Lee V W. Scattering of anti-plane SH-waves by a shallow semi-elliptical hill with a concentric elliptical tunnel. Journal of Earthquake Engineering, 2016, 20(3): 363–382

[17] Kara H F. Diffraction of plane SH-wave by a cylindrical cavity in an infinite wedge. Procedia Engineering, 2016, 161: 1601–1607

[18] Gao Y, Dai D, Zhang N, Wu Y, Mahfouz A H. Scattering of plane and cylindrical SH-waves by a horseshoe shaped cavity. Journal of Earthquake and Tsunami, 2016, 10(2): 1–23

[19] Wang G, Liu D. Scattering of SH-wave by multiple circular cavities in half space. Earthquake Engineering and Engineering Vibration, 2002, 1(1): 36–44

[20] Shi W P, Liu D K, Song Y T, Chu J L, Hu A Q. Scattering of circular cavity in right-angle planar space to steady SH-wave. Applied Mathematics and Mechanics, 2006, 27(12): 1619–1626

[21] Lee V W, Manoogian M E, Chen S. Anti-plane SH-deformations near a surface rigid foundation above a subsurface rigid circular tunnel. Earthquake Engineering and Engineering Vibration, 2002, 1(1): 27–35

[22] Smerzini C, Aviles J, Paolucci R, Sánchez-Sesma F J. Effect of underground cavities on surface earthquake ground motion under SH-wave propagation. Earthquake Engineering & Structural Dynamics, 2009, 38(12): 1441–1460

[23] Liu D K, Lin H. Scattering of SH-waves by an interacting interface linear crack and a circular cavity near bimaterial interface. Chinese Journal of Theoritical and Applied Mechanics, 2004, 20(3): 317–326

[24] Tsaur D H, Chang K H. Multiple scattering of SH-waves by an embedded truncated circular cavity. Journal of Marine Science and Technology, 2012, 20(1): 73–81

[25] He C, Zhou S, Guo P, Di H, Zhang X. Analytical model for vibration prediction of two parallel tunnels in a full-space. Journal of Sound and Vibration, 2018, 423: 306–321

[26] Lee K M, Rowe R K. An analysis of three-dimensional ground movements: The Thunder Bay tunnel. Canadian Geotechnical Journal, 1991, 28(1): 25–41

[27] Molinero J, Samper J, Juanes R. Numerical modeling of the transient hydrogeological response produced by tunnel construction in fractured bedrocks. Engineering Geology, 2002, 64(4): 369–386

[28] Karakus M, Ozsan A, Başarır H. Finite element analysis for the twin metro tunnel constructed in Ankara Clay, Turkey. Bulletin of Engineering Geology and the Environment, 2007, 66(1): 71–79

[29] Yamamoto K, Lyamin A V, Wilson D W, Sloan S W, Abbo A J. Stability of dual circular tunnels in cohesive-frictional soil subjected to surcharge loading. Computers and Geotechnics, 2013, 50: 41–54

[30] Zhang Z X, Liu C, Huang X, Kwok C Y, Teng L. Three-dimensional finite-element analysis on ground responses during twin-tunnel construction using the URUP method. Tunnelling and Underground Space Technology, 2016, 58: 133–146

[31] Das R, Singh P K, Kainthola A, Panthee S, Singh T N. Numerical analysis of surface subsidence in asymmetric parallel highway tunnels. Journal of Rock Mechanics and Geotechnical Engineering, 2017, 9(1): 170–179

[32] Tsinidis G. Response of urban single and twin circular tunnels subjected to transversal ground seismic shaking. Tunnelling and Underground Space Technology, 2018, 76: 177–193

[33] Yiouta-Mitra P, Kouretzis G, Bouckovalas G, Sofianos A. Effect of underground structures in earthquake resistant design of surface structures. Dynamic Response and Soil Properties, 2007, 1–10

[34] Chakeri H, Hasanpour R, Hindistan M A, Unver B. Analysis of interaction between tunnels in soft ground by 3-D numerical modeling. Bulletin of Engineering Geology and the Environment, 2011, 70(3): 439–448

[35] Besharat V, Davoodi M, Jafari M K. Effect of underground structures on free-field ground motion during earthquakes. In: The 15th World Conference on Earthquake Engineering. Lisbon: SPES, 2012

[36] Do N A, Dias D, Oreste P, Djeran-Maigre I. Three-dimensional numerical simulation of a mechanized twin tunnels in soft ground. Tunnelling and Underground Space Technology, 2014, 42: 40–51

[37] Jamshasb F, Shahrbabaki M K. The impact of approach in shallow twin tunnels on the convergence. Procedia Earth and Planetary Science, 2015, 15: 325–329

[38] Narayan J P, Kumar D, Sahar D. Effects of complex interaction of Rayleigh-waves with tunnel on the free surface ground motion and the strain across the tunnel-lining. Natural Hazards, 2015, 79(1): 479–495

[39] Ziaei A, Ahangari K. The effect of topography on stability of shallow tunnels case study: The diversion and conveyance tunnels of Safa Dam. Transportation Geotechnics, 2018, 14: 126–135

[40] Beskos D E. Boundary element methods in dynamic analysis. Applied Mechanics Reviews, 1987, 40(1): 1–23

[41] Luco J E, de Barros F C P. Dynamic displacements and stresses in the vicinity of a cylindrical cavity embedded in a half-space. Earthquake Engineering & Structural Dynamics, 1994, 23(3): 321–340

[42] Yu M C, Dravinski M. Scattering of a plane harmonic SH-wave by a completely embedded corrugated scatterer. International Journal for Numerical Methods in Engineering, 2009, 78(2): 196–214

[43] Parvanova S L, Dineva P S, Manolis G D, Wuttke F. Seismic response of lined tunnels in the half-plane with surface topography. Bulletin of Earthquake Engineering, 2014, 12(2): 981–1005

[44] Liu Z, Liu L. An IBEM solution to the scattering of plane SH-waves by a lined tunnel in elastic wedge space. Earth Science, 2015, 28(1): 71–86

[45] Benites R, Aki K, Yomogida K.Multiple scattering of SH-waves in 2-D media with many cavities. Pure and Applied Geophysics, 1992, 138(3): 353–390

[46] Panji M, Asgari Marnani J, Tavousi Tafreshi S. Evaluation of effective parameters on the underground tunnel stability using BEM. Journal of Structural Engineering and Geotechnic, 2011, 1(2): 29–37

[47] Panji M, Koohsari H, Adampira M, Alielahi H, Asgari Marnani J. Stability analysis of shallow tunnels subjected to eccentric loads by a boundary element method. Journal of Rock Mechanics and Geotechnical Engineering, 2016, 8(4): 480–488

[48] Panji M, Ansari B. Modeling pressure pipe embedded in two-layer soil by a half-plane BEM. Computers and Geotechnics, 2017, 81(C): 360–367

[49] Takemiya H, Fujiwara A. SH-wave scattering and propagation analyses at irregular sites by time domain BEM. Bulletin of the Seismological Society of America, 1994, 84(5): 1443–1455

[50] Kamalian M, Jafari M K, Ghayamghamian M R, Shafiee A, Hamzehloo H, Haghshenas E, Sohrabi-Bidar A. Site effect microzonation of Qom, Iran. Engineering Geology, 2008, 97(1–2): 63–79

[51] Rice J M, Sadd M H. A note on computing elastodynamic full field displacements arising from subsurface singular sources. Mechanics Research Communications, 1984, 11(6): 385–390

[52] Belytschko T, Chang H S. Simplified direct time integration boundary element method. Journal of Engineering Mechanics, 1988, 114(1): 117–134

[53] Hirai H. Analysis of transient response of SH-wave scattering in a half space by the boundary element method. Engineering Analysis, 1988, 5(4): 189–194

[54] Panji M, Kamalian M, Asgari-Marnani J, Jafari M K. Transient analysis of wave propagation problems by half-plane BEM. Geophysical Journal International, 2013, 194(3): 1849–1865

[55] Panji M, Kamalian M, Marnani J A, Jafari M K. Anti-plane seismic response from semi-sine shaped valley above embedded truncated circular cavity: A time-domain half-plane BEM. International Journal of Civil Engineering, 2014, 12(2): 193–206

[56] Panji M, Kamalian M, Marnani J A, Jafari M K. Analyzing seismic convex topographies by a half-plane time-domain BEM. Geophysical Journal International, 2014, 197(1): 591–607

[57] Panji M, Ansari B. Transient SH-wave scattering by the lined tunnels embedded in an elastic half-plane. Engineering Analysis with Boundary Elements, 2017, 84: 220–230

[58] Panji M, Mojtabazadeh-Hasanlouei S. Time-history responses on the surface by regularly distributed enormous embedded cavities: Incident SH-waves. Earth Science, 2018, 31: 1–17

[59] Panji M, Mojtabazadeh-Hasanlouei S. Seismic amplification pattern of the ground surface in presence of twin unlined circular tunnels subjected to SH-waves. Journal of Transportation Infrastructure Engineering, 2019, 5(3): 111–134 (in Persian)

[60] Huang L, Liu Z, Wu C, Liang J. The scattering of plane P, SV-waves by twin lining tunnels with imperfect interfaces embedded in an elastic half-space. Tunnelling and Underground Space Technology, 2019, 85: 319–330

[61] Panji M, Mojtabazadeh-Hasanlouei S, Yasemi F. A half-plane time-domain BEM for SH-wave scattering by a subsurface inclusion. Computers & Geosciences, 2020, 134: 104342

[62] Mojtabazadeh-Hasanlouei S, Panji M, Kamalian M. On subsurface multiple inclusions model under transient SH-wave propagation. Waves in Random and Complex Media, 2020

[63] Panji M, Mojtabazadeh-Hasanlouei S. Transient response of irregular surface by periodically distributed semi-sine shaped valleys: Incident SH-waves. Journal of Earthquake and Tsunami, 2020, 14(1): 2050005

[64] Chen L, Liu C, Zhao W, Liu L. An isogeometric approach of two-dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution. Computer Methods in Applied Mechanics and Engineering, 2018, 336: 507–532

[65] Chen L, Lu C, Lian H, Liu Z, Zhao W, Li S, Chen H, Bordas S P A. Acoustic topology optimization of sound absorbing materials directly from subdivision surfaces with isogeometric boundary element methods. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112806

[66] Shaaban A M, Anitescu C, Atroshchenko E, Rabczuk T. Shape optimization by conventional and extended isogeometric boundary element method with PSO for two-dimensional Helmholtz acoustic problems. Engineering Analysis with Boundary Elements, 2020, 113: 156–169

[67] Shaaban A M, Anitescu C, Atroshchenko E, Rabczuk T. Isogeometric boundary element analysis and shape optimization by PSO for 3D axi-symmetric high frequency Helmholtz acoustic problems. Journal of Sound and Vibration, 2020, 486: 115598

[68] Kolymbas D. Tunneling and Tunnel Mechanics. Berlin: Springer, 2005

[69] Morse P M, Feshbach H. Methods of Theoretical Physics. New York: McGraw-Hill Book Company, 1953

[70] Eringen A C, Suhubi E S. Elastodynamics. New York: Academic Press, 1975

[71] Brebbia C A, Dominguez J. Boundary Elements: An Introductory Course. Southampton: Computational Mechanics Publications, 1989

[72] Dominguez J. Boundary Elements in Dynamics. Southampton: Computational Mechanics Publications, 1993

[73] Panji M, Ansari B. Anti-plane seismic ground motion above twin horseshoe-shaped lined tunnels. Innovative Infrastructure Solutions, 2020, 5(1): 1–16

[74] Ji D, Lei W, Li H. Corner treatment by assigning dual tractions to every node for elastodynamics in TD-BEM. Applied Mathematics and Computation, 2016, 284: 125–135

[75] Mathworks. MATLAB, The Language of Technical Computing V.9.9. Massachusetts: The MathWorks Inc, 2000

[76] Trifunac M D. Scattering of plane SH-waves by a semi-cylindrical canyon. Earthquake Engineering & Structural Dynamics, 1972, 1(3): 267–281

[77] Borcherdt R D. Estimates of site-dependent response spectra for design (methodology and justification). Earthquake Spectra, 1994, 10(4): 617–653

RIGHTS & PERMISSIONS

Higher Education Press