Transient Dynamic Response and Anti-Seismic Measures of Deep-Buried Composite Lining Tunnels Subject to Blasting SV-Wave Disturbance

1  Introduction

Drilling and blasting, as an efficient and economical rock excavation technique, is widely used in underground chamber excavation, mining engineering, and other fields [1]. Nowadays, a growing number of newly constructed tunnels are situated in close proximity to existing ones, such as during the expansion and renovation of railway tunnels or in the construction of intersecting tunnels in metro projects [2]. These blasting activities may pose threats to the structural safety of existing tunnels, leading to issues like lining cracks, deformation, or even collapse [3,4]. Therefore, an in-depth investigation into the dynamic response of tunnel linings to blast-induced seismic waves is critical for accurately evaluating tunnel safety and optimizing engineering efficiency.

Lining failure is primarily attributed to differential displacements at various locations caused by stress [5]. Hence, stress and displacement will be the key indicators examined in this study. Pao et al. [6] systematically investigated the dynamic stress concentration around a single cavity and a circular lined tunnel in an infinite space under incident elastic waves. Subsequent researchers further found that under certain frequencies of incident waves, the maximum hoop stress concentration induced by incident SV-waves can be 40%–60% higher than that caused by P-waves [7]. Lu et al. [8] simplified blast-induced seismic waves as plane harmonic waves to study the scattering problem around tunnels and examined the effects of the lining’s shear modulus, wave number, and thickness on the dynamic stress concentration factor (DSCF). Li et al. [9] employed the wave function expansion method to derive expressions for the DSCF in surrounding rock and linings under blast-induced seismic waves in a full space, focusing primarily on the steady-state response (i.e., frequency domain). Later, researchers introduced the large-arc assumption to extend these solutions to half-space problems and derived series solutions for the scattering of P-, SV-, and SH-waves by a single cavity or a circular lined tunnel in a half-space [10–12]. In recent years, Ding and Tong investigated the scattering and stress concentration phenomena of plane waves acting on a circular lining in a half-space based on nonlocal Biot theory [13]. Subsequently, they applied auxetic materials with negative Poisson’s ratio [14] and locally resonant metamaterials [15] to seismic vibration isolation through numerical and experimental approaches. Furthermore, the vibration isolation characteristics of nonlinear metasurfaces in saturated porous media were analyzed theoretically [16].

In practical tunnel engineering, due to factors such as waterproofing membranes or cracks, the surrounding rock and lining are not perfectly bonded. Thus, intense blasting or seismic events may easily cause fracture or collapse of the lining at these weakly bonded areas. Such weak interfaces can be regarded as imperfect [17]. In recent years, some researchers have introduced linear spring models [18,19] to describe imperfect interfaces, assuming stress continuity but displacement discontinuity. Fang et al. [20] studied the effect of an imperfect interface on the stress concentration in the lining of a twin-tunnel under blast P-waves. Yi et al. [19] used the wave function expansion method with a spring model to simulate an imperfect interface, showing that changes in interface stiffness significantly affect the distribution of the DSCF in the surrounding rock and lining. More recently, Xiang et al. [21] introduced an elastic-slip interface to characterize imperfect bonding.

It is important to note that the aforementioned studies are limited to steady-state incident waves. In reality, engineering activities such as blasting and drilling generate transient waves, which propagate through the surrounding rock and can cause damage to both the rock mass and existing structures [22]. Therefore, building on steady-state response analyses, some scholars have simplified blast loads as triangular, half-sine, or exponential functions, and used Fourier or Laplace transforms to study the transient response of circular tunnels under blast P-waves [23]. While a half-sine wave load features equal rise and decay times, and an exponential model is often mathematically complex and less tractable, the triangular load model offers a simpler form with independent control over key parameters such as peak pressure, rise time, and duration. These parameters have clear physical interpretations, allowing the model to effectively capture the essential characteristics of a blast load, making it a widely adopted idealization [23,24]. In recent years, researchers have increasingly turned to the time-domain boundary element method (TD-BEM) and the dual reciprocity boundary element method (DR-BEM) to model transient incident waves [25–27]. Existing research on transient response has mainly focused on the dynamic behavior of cavities and single-layer linings under blast P-waves, largely overlooking the significant impact of blast SV-waves on underground structures. While these methods significantly reduce computational time, they are unable to accurately capture the distinctive characteristics of blasting loads. In light of this, this study focuses on composite linings, which are widely used under complex geological conditions and strong dynamic loads. The interaction between their multiple structural layers leads to a more complex dynamic response mechanism compared to single-layer linings, underscoring the specificity and importance of this research.

In summary, this study aims to investigate the transient dynamic characteristics of composite tunnels and explore effective seismic mitigation measures through an analysis under blast-induced SV-wave loading. First, analytical solutions for the transient stress and displacement responses in the surrounding rock and lining are derived using the wave function expansion method and Fourier transform. Then, the effects of blast SV-waves on the transient response of the surrounding rock and lining are compared through numerical calculations in MATLAB. Finally, the influence of parameters such as blast loading duration, lining thickness, and elastic modulus on the transient response of the tunnel lining is analyzed.

2  Computational Model and Steady-State Response Solution

2.1 Computational Model

The shock wave generated by the explosion of explosives in rock attenuates into a stress wave at a distance of 125–150 times the charge radius (r) from the source, and further decays into a seismic wave beyond 150r [8]. In the far-field region, the blast-induced seismic wave can be approximated as a plane wave. Taking the plane SV-wave as an example, when it propagates perpendicularly to the tunnel axis, the problem can be simplified as a plane strain case. The computational model of the composite lining tunnel, the characteristics of the blast load, and the corresponding coordinate system are shown in Fig. 1. Here, the black dashed line indicates the division between the near-field and mid-field blast zones, while the blue dashed line represents the boundary between the mid-field and far-field blast zones. The area beyond the blue line—the far-field blast zone—is the primary focus of this study. The medium is assumed to be an isotropic linear elastic material. In the figure, (x, y) denotes the Cartesian coordinate system, while (θ, r) represents the polar coordinate system; r2 and r1 indicate the inner and outer radii of the initial support, respectively, and r3 represents the inner radius of the secondary lining.

Figure 1: Schematic diagram of the calculation model and blast load characterization

The time-domain characteristics of blast-induced seismic waves are highly complex during the blasting process. Therefore, to simplify the theoretical analysis, the blast-induced seismic wave is idealized as a triangular pulse, as illustrated in Fig. 1, where f(t) denotes the time history of the blast load, fm represents the peak value of the seismic wave, and tr and ts correspond to the rise time and total duration of the blast load, respectively. The transient response from any transient wave can be determined by the superposition of its harmonic components. Therefore, the theoretical formulation for the dynamic response of a composite lined tunnel under harmonic SV-wave excitation is derived first.

2.2 Wave Field Analysis

According to the wave function expansion method, the incident wave field of a plane SV-wave propagating in the positive x-direction can be expressed as [6]:

ψ(i)=φ0∑n=0∞εninJn(β1r)sin⁡(nθ)e−iωt(1)

where: φ0 denotes the amplitude of the incident wave, β1 is the wavenumber of the SV-wave, ω represents the angular frequency of the incident wave, Jn refers to the Bessel function of the first kind of order n. When n = 0 occurs, εn=1 follows; when n≥1 occurs, εn=2 results.

When an SV-wave impinges on the interface between the surrounding rock and the initial support, it is reflected as both P- and SV-waves. The corresponding potential functions for the P-wave and the SV-wave can be expressed respectively as:

{φ(r)=∑n=0∞AnHn(1)(α1r)cos⁡(nθ)e−iωtψ(r)=∑n=0∞BnHn(1)(β1r)sin⁡(nθ)e−iωt(2)

Therefore, the total wave field potential in the surrounding rock can be expressed as the superposition of the incident wave potential and all scattered wave potentials:

{φ1=φ(r)ψ1=ψ(i)+ψ(r)(3)

Similarly, the total wave fields within the initial support and the secondary lining can be expressed respectively as:

{φ2=∑n=0∞[CnHn(1)(α2r)+DnHn(2)(α2r)]cos⁡(nθ)ψ2=∑n=0∞[FnHn(1)(β2r)+GnHn(2)(β2r)]sin⁡(nθ)φ3=∑n=0∞[KnHn(1)(α3r)+LnHn(2)(α3r)]cos⁡(nθ)ψ3=∑n=0∞[MnHn(1)(β3r)+NnHn(2)(β3r)]sin⁡(nθ)(4)

where: αm,βm(m=1,2,3) represent the wavenumbers of P-waves and SV-waves in the initial support and the secondary lining, respectively; Hn(1) and Hn(2) denote the Hankel functions of the first and second kinds of order n. An,Bn,Cn,Dn,Fn,Gn,Kn,Ln,Mn,Nn are undetermined coefficients.

2.3 Solution of Steady-State Response

When a harmonic SV-wave is incident upon the interfaces between the surrounding rock, initial support, and secondary lining, the resulting reflected and transmitted waves give rise to specific displacement and stress fields. The corresponding expressions for these fields are derived and presented below:

{ur=∂φ∂r+1r∂ψ∂θuθ=1r∂φ∂θ−∂ψ∂rσrr=λ∇2φ+2μ[∂2φ∂r2+∂∂r(1r∂ψ∂θ)]σθθ=λ∇2φ+2μ[1r(∂φ∂r+1r∂2φ∂θ2)+1r(1r∂ψ∂θ−∂∂r∂θ)]σrθ=2μ(1r∂2φ∂r∂θ−1r2∂φ∂θ)+μ[1r2∂2ψ∂θ2−r∂∂r(1r∂ψ∂r)](5)

where: ∇2φ=∂2φ∂r2+1r∂φ∂r+1r2∂2φ∂θ2; ∇2 denotes the Laplace operator, ur,uθ represent the radial and tangential displacements in polar coordinates, respectively;

σrr,σθθ and σrθ correspond to the radial, tangential, and shear dynamic stresses, respectively; λ,μ represents the Lamé constants of the medium.

Based on the displacement continuity condition at the interface between the surrounding rock and the initial support, the following relation holds at r=r1:

{σrr(2)=σrr(1),σrθ(2)=σrθ(1)ur(2)=ur(1),uθ(2)=uθ(1)(6)

At the contact surface r=r2 between the initial support and the secondary lining, the following condition is satisfied:

{σrr(3)=σrr(2),σrθ(3)=σrθ(2)ur(3)=ur(2),uθ(3)=uθ(2)(7)

At the inner surface of the tunnel (r=r3), the stress-free boundary condition applies:

σrr(3)=0,σrθ(3)=0(8)

Substituting Eqs. (3) and (4) into Eq. (5), and then incorporating the resulting expressions for displacement and stress into Eqs. (6) to (8), the following matrix Eq is obtained:

{Aij}[Xi]=[Mj],i,j=1,2,…10(9)

where: X1=An,X2=Bn⋅⋅⋅X9=Mn,X10=Nn. By applying the orthogonality of trigonometric functions, the coefficients {Aij} and [Mj] can be determined. The detailed expressions for their individual elements are provided in Appendix A. By solving the matrix Eq, the undetermined coefficients in the potential functions of each medium can be obtained, thereby yielding the expression for the DSCF in each domain:

DSCF=σθθ∗=|σθθ/σ0|(10)

where: σ0=−μβ12φ0 representing the maximum stress value of the incident wave, μ represents the Lamé constant.

3  Solution for Transient Response

The previous section analyzed the steady-state response of the composite lined tunnel under harmonic SV-wave excitation. However, in engineering practice, blasting exerts transient dynamic disturbances on underground structures. By employing Fourier transform techniques, the blast-induced waveform can be decomposed into its harmonic components. The transient response expression can then be derived through superposition of these components:

g(xi,t)=12π∫−∞+∞χ(x,ω)F(ω)e−iωtdω(11)

where: χ(xi,ω) represents the admittance function derived from the steady-state response in Eq. (10); F(ω) denotes the Fourier transform of the blast load time-history.

Referring to existing literature [28], the blast load is defined by a trigonometric function, as given in Eq. (12):

f(t)={0t<0t/tt0≤t<tr(ts−t)/(ts−tr)tr≤t<ts0t≥ts(12)

In the computational model shown in Fig. 1, the time origin is set when the incident wave first reaches the outer boundary of the initial support (i.e., at r = r1). The elapsed time t is normalized by the time required for a wave to travel across the outer radius of the initial support, as given by Eq. (13):

τ=cpt/r1(13)

By introducing the δ-function and the Heaviside step function, the integration process in Eq. (13) can be simplified. The Fourier transform of the blast load, δ(t), is given by:

12π∫−∞+∞δ(t)e−iωtdt=12π(14)

Therefore, the unit impulse response induced by the δ(t)-function can be expressed as:

gδ(x,t)=12π∫−∞+∞χ(x,ω)e−iωtdω(15)

when the blast load propagates to the boundary of the initial support, the integration of Eq. (15) yields:

gh(x,t)=∫0tgδ(x,t)dt(16)

Substituting the expanded time factor from Eq. (15) into Eq. (16) yields the expression for the Heaviside step function:

gh(x,t)=2π∫0∞χ(x,ω)ωsin⁡(ωt)dω(17)

The transient response for an arbitrary function can be expressed based on the principle of Duhamel’s integral as:

g(x,t)=f(0)gh(t)+∫0tf′(τ)gh(t−τ)(18)

Substituting Eqs. (12) and (17) into Eq. (18) yields the expression for the transient response around the composite lined tunnel.

When 0≤t<tr occurs:

g(xi,t)=∫0t1trdτ2π∫0∞R(ω)sin⁡ω(t−τ)ωdω=2πtr∫0∞R(ω)(1−cos⁡ωt)ω2dω(19a)

when tr≤t<ts occurs:

g(xi,t)=2πtr∫0∞R(ω)[cos⁡ω(t−tr)−cos⁡ωt]ω2dω−2π(ts−tr)∫0∞R(ω)[1−cos⁡ω(t−tr)]ω2dω(19b)

when t≥ts occurs:

g(xi,t)=2πtr∫0∞R(ω)[cos⁡ω(t−tr)−cos⁡ωt]ω2dω−2π(ts−tr)∫0∞R(ω)[cos⁡ω(t−ts)−cos⁡ω(t−tr)]ω2dω(19c)

Theoretically, the steady-state response of the composite lined tunnel can be calculated using Eq. (9). However, obtaining the analytical expression for the admittance function χ(x,ω) is challenging, which complicates the direct integration of Eq. (19). Therefore, this study employs the trapezoidal quadrature method [21] for numerical approximation. In this method, the real part of the steady-state response for each wavenumber is first determined by Eq. (9). The transient response contribution from each segment, as expressed in Eq. (19), is then calculated employing the trapezoidal rule. The overall transient response of the system is finally obtained by superimposing frequency components through Duhamel integration, which is numerically implemented in the computational software using the trapz function.

4  Results and Analysis

4.1 Verification of Results

To evaluate the reliability of the theoretical model, the physical parameters in this study were first set to be consistent with those in existing literature and compared with the results reported by Fan et al. [5]. For ease of comparison with previous studies, the following dimensionless parameters were adopted for the case of harmonic SV-wave incidence: μ~=2.9 (ratio of shear moduli between the surrounding rock and the lining), γ=1.5 (ratio of P-wave velocities in the surrounding rock and the lining), v1=0.25 (Poisson’s ratio of the surrounding rock), v2=0.2 (Poisson’s ratio of the lining, for single-layer case), η=1.2 (ratio of outer to inner radius of the lining, for single-layer case), and β1r3 (dimensionless wavenumber). The computational results are shown in Fig. 2. The degenerate solution derived in this study agrees well with the numerical results from Ref. [5], confirming that the proposed theoretical model can be reliably applied to further analyses.

Figure 2: Comparison of the degradation solution in this article with the results in the Ref. [5]

4.2 Case Study and Analysis

The blasting construction for the Hongshan South Road tunnel group (delineated by the blue dashed box in Fig. 3) poses potential impacts to the operational safety of the existing Metro Line 1 tunnel (indicated by the red dashed box). Therefore, this tunnel group is selected as a case study for parametric analysis in this paper. The relative position and distance between them are shown in Fig. 3 below. Referring to existing Ref. [9], the computational parameters in this study are selected as follows: the inner radius of the secondary lining r3 is set to 5 m, while other geometric and material parameters of the tunnel are summarized in Table 1.

Figure 3: Relative position of tunnels in Hongshan South Road

4.2.1 Distribution Patterns of Transient Dynamic Response in Composite Lined Tunnels

Fig. 4 shows the circumferential distribution of the peak dynamic stress concentration factor (DSCFmax) at the inner and outer boundaries of the tunnel lining. The ratio ts/tr characterizes the unloading rate under blast loading. A smaller value of this ratio indicates a faster unloading process; when ts/tr = 1 it represents instantaneous unloading of the blast load. For the blast loading parameters in Fig. 4, ts/tr = 5, tr = 5. As clearly observed in the figure, the distribution patterns of DSCFmax at the inner and outer boundaries are generally consistent during blast loading, exhibiting a “four-leaf clover” shape. Significant dynamic stress concentration occurs at θ = 45°, 135°, 225° and 315°. Moreover, the DSCFmax at the inner boundary of the composite lining is greater than that at the outer boundary. This occurs because the higher elastic modulus of the secondary lining results in a greater wave impedance, leading to the reflection and superposition of stress waves as they transmit from the primary lining, which generates a dynamic amplification effect. Furthermore, the inner boundary of the secondary lining is a free surface, lacking external confinement, resulting in a more critical stress state.

Figure 4: Distribution curves of DSCFmax at the inner and outer boundaries of the composite lining

4.2.2 Influence of Blast Loading Duration on the Dynamic Stress Concentration Factor

Fig. 5 shows the time-history of DSCF under different blast loading durations. The computational parameters align with previous sections. A triangular blast load is defined by parameters ts/tr = 5 and tr. Fig. 5 presents the DSCF time histories in the secondary lining and surrounding rock under SV-wave excitation for varying blast loading durations. Generally, the transient DSCF rapidly rises to a positive peak, followed by a sharp decline to negative values before gradually attenuating to zero. This behavior reflects the transition from compressive to tensile stress states during wave scattering around the lining. For the monitoring point at θ = 90°, the blast loading duration has a noticeable influence on the DSCF in both the secondary lining and the surrounding rock, although the peak values remain relatively similar. In contrast, at the θ = 45° location, the peak DSCF in the secondary lining is approximately twice that observed in the surrounding rock, a trend consistent with previous findings. Furthermore, the influence of blast loading duration at this location is considerably less pronounced. The DSCF peak exhibits a negative correlation with the blast loading duration under transient SV-wave action. Shorter durations result in higher DSCF amplitudes due to the more rapid input and release of energy, which consequently induces a more pronounced dynamic stress concentration and structural response.

Figure 5: Effect of blast loading duration on transient DSCF: (a) θ = 45° rock; (b) θ = 45° secondary lining; (c) θ = 90° rock; (d) θ = 90° secondary lining

4.3 Optimization of Lining Parameters

The circumferential distribution and time-history characteristics indicate that the peak DSCF is predominantly concentrated at the haunch regions in both the surrounding rock and the secondary lining. Therefore, the analysis in this section focuses on DSCFmax at point θ = 45°. The parametric analysis of this paper will focus on the elastic modulus and thickness, two parameters that play a leading role in the dynamic response of the lining and are most critical in the engineering design. By analyzing the physico-mechanical parameters of the Hongshan South Road Tunnel, and employing the method of determining undetermined coefficients, this study investigates the variation of peak DSCF in the surrounding rock and the inner lining under different values of shear modulus and lining thickness. Such an analysis aids in identifying the most influential parameters on the dynamic response of the tunnel structure, thereby providing guidance for tunnel design and construction. The blast loading parameters are set as ts/tr = 5 and tr = 5, while the remaining computational parameters are consistent with those previously defined.

4.3.1 Elastic Modulus of the Lining

As illustrated in Fig. 6, under transient SV-wave loading, the elastic moduli of the initial support and secondary lining exhibit generally consistent influences on the DSCFmax of both the surrounding rock and the secondary lining. Specifically, the DSCFmax of the secondary lining shows a positive correlation with both the elastic modulus of the initial support (E2) and its own elastic modulus (E3), while the DSCFmax of the surrounding rock is negatively correlated with E2 and E3. Fig. 6a presents the variation in DSCFmax at point θ = 45° for the secondary lining and surrounding rock as E2 changes, while E3 is held constant at 32 GPa. Overall, Fig. 6a indicates that the DSCFmax of the secondary lining exhibits a positive correlation with the elastic modulus of both the primary lining (E2) and itself (E3), whereas the DSCFmax of the surrounding rock shows a negative correlation with E2 and E3. When E2 is small, the DSCFmax of the secondary lining shows significant variation. However, when E2 is large, the DSCFmax exhibits noticeable oscillations. This is likely due to the higher shear wave velocity in the high-modulus primary lining, which causes the SV-wave energy to concentrate and local stresses to change rapidly, thereby inducing strong oscillatory effects. Fig. 6b displays the variation in DSCFmax at point θ = 45° with E3, while E2 remains constant at 25 GPa. The results indicate that the DSCFmax of the secondary lining increases significantly with E3—rising from 2.0 to 3.9—while the DSCFmax of the surrounding rock decreases with E3. The change in the DSCFmax of the secondary lining is notably more sensitive to variations in E3. By comparing Fig. 6a,b, it can be concluded that the dynamic response on both the inner and outer surfaces of the lining is more significantly affected by the elastic modulus of the secondary lining. Therefore, in seismic mitigation design, while ensuring the stability of the initial support, selecting a material with a lower elastic modulus for the secondary lining can effectively reduce dynamic stress concentration.

Figure 6: Effect of elastic modulus on DSCFmax: (a) Elastic modulus of initial support/GPa; (b) Elastic modulus of secondary lining/GPa

4.3.2 Lining Thickness

As shown in Fig. 7, under transient SV-wave excitation, the thicknesses of both the initial support and the secondary lining exhibit a negative correlation with the DSCFmax in the surrounding rock and the secondary lining. Fig. 7a shows the variation of DSCFmax in the secondary lining and surrounding rock at point θ = 90° with the primary lining thickness, while the secondary lining thickness remains constant at 0.45 m. It can be observed that increasing the primary lining thickness significantly improves the internal stress state of the lining. The DSCFmax of the surrounding rock is more sensitive to changes in the primary lining thickness, decreasing from 1.58 to 1.04, a reduction of approximately 34%. Fig. 7b shows the variation of DSCFmax at point θ = 45° with the thickness of the secondary lining, while the initial support thickness remains constant at 0.15 m. It can be observed that the DSCFmax of both the surrounding rock and the secondary lining decreases as the thickness of the secondary lining increases. Notably, the DSCFmax of the surrounding rock shows a substantial reduction—from 2.00 to 1.07, representing a decrease of approximately 50%—demonstrating a high sensitivity to changes in the secondary lining thickness. These findings indicate that the thickness of the secondary lining has a considerable influence on the DSCFmax on both the inner and outer sides of the lining. Therefore, in seismic mitigation design, the thickness of the secondary lining should be carefully considered to effectively control dynamic stress concentration.

Figure 7: Effect of lining thickness on DSCFmax: (a) Initial support thickness/m; (b) Secondary lining thickness/m

5  Conclusion

Based on the wave function expansion method and Fourier transform, this study derived an analytical solution for the transient response of a circular composite-lined tunnel under incident transient SV-waves. The computational procedure was validated by comparing the degenerate steady-state solution with existing research. Following this validation, and by integrating practical engineering conditions with lining parameter optimization, the following conclusions were drawn:

1.   In tunnel seismic design, special attention must be paid to the haunch area on the inner side of the secondary lining, as a short rise time of blasting loads will cause a significant increase in the peak stress in this region, necessitating enhanced protective measures.

2.   An increase in the stiffness of the primary support raises the DSCF at the surrounding rock boundary while reducing it at the inner boundary of the secondary lining; conversely, increasing the stiffness of the inner boundary of the secondary lining has the opposite effect. Therefore, on the premise of ensuring the stability of the primary support, secondary lining materials with lower stiffness should be selected in the design.

3.   Increasing the lining thickness of either the primary support or the secondary lining can substantially reduce the DSCF within the lining, a principle that should be incorporated into the seismic resistance and mitigation design of composite linings.

Future research could consider more complex tunnel geometries, nonlinear constitutive models of the surrounding rock, and three-dimensional spatial effects to more accurately simulate real-world engineering scenarios.

Acknowledgement: The authors gratefully acknowledge the dedication and contributions of all team members.

Funding Statement: The work was supported by the Research Project on Micro-Vibration Blasting Technology for Tunnels in High-Altitude Cold Regions (2024HX01), the Jiangxi “Ganpo Jun Cai” Program for Young Sci-Tech Talents (2024QT04) and the Natural Science Foundation of Jiangxi Province (20242BAB204095).

Author Contributions: The individual contributions of the authors are specified as follows: Qunjie Huang: Writing—original draft, Methodology, Investigation, Funding acquisition. Yu Huang: Supervision, Funding acquisition, Conceptualization. Yangqing Liu: Supervision, Investigation. Qiaoming Guo: Writing—review & editing, Methodology. Zhiyun Liu: Methodology, Validation, Resources, Data curation. Haibin Ding: Supervision, Funding acquisition. Lihua Li made critical revisions to the manuscript. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

Appendix A

ε11(l)(αr)=(n2+n−1/2β2r2)Qn(l)(αr)−αrQn−1(l)(αr)ε12(l)(βr)=n[−(n+1)Qn(l)(βr)+βrQn−1(l)(βr)]ε21(l)(αr)=−(n2+n+1/2β2r2−α2r2)Qn(l)(αr)+αrQn−1(l)(αr)ε22(l)(βr)=n[(n+1)Qn(l)(βr)−βrQn−1(l)(βr)]ε41(l)(αr)=−n[−(n+1)Qn(l)(αr)+αrQn−1(l)(αr)]ε42(l)(βr)=−(−n2+n−1/2β2r2)Qn(l)(βr)+βrQn−1(l)(βr)ε71(l)(αr)=αrQn−1(l)(αr)−nQn(l)(αr)ε72(l)(βr)=nQn(l)(βr)ε81(l)(αr)=−nQn(l)(αr)ε82(l)(βr)=−[βrQn−1(l)(βr)−nQn(l)(βr)]Qn(l)(⋅)={Jn(⋅),l=1Hn(1)(⋅),l=3Hn(2)(⋅),l=4

{M1=inφ0εnE411(β1r)M2=inφ0εnE111(β1r)M3=inφ0εnE811(β1r)M4=inφ0εnE711(β1r)

{A11=−E113(α1r),A12=−E123(β1r)A13=μ2E123(α2r)A14=μ2E124(α2r)A15=E123(β2r),A16=μ2E124(β2r),A17=A18=A19=A110=0A21=−E413(α1r),A22=−E423(β1r)A23=μ2E423(α2r)A24=μ2E424(α2r)A25=E423(β2r),A26=μ2E424(β2r),A27=A28=A29=A210=0A31=−E713(α1r),A32=−E723(β1r)A33=μ2E723(α2r)A34=μ2E724(α2r)A35=E723(β2r),A36=μ2E724(β2r),A37=A38=A39=A310=0A41=−E813(α1r),A42=−E823(β1r)A43=μ2E823(α2r)A44=μ2E824(α2r)A45=E823(β2r),A46=μ2E824(β2r),A47=A48=A49=A410=0

{A51=A52=0,A53=μ2E123(α2r),A54=μ2E124(α2r),A55=E123(β2r),A56=μ2E124(β2r)A57=−μ1E113(α3r),A58=−μ1E114(α3r),A59=−μ1E123(β3r),A510=−μ1E124(β3r)A61=A62=0,A63=μ2E423(α2r),A64=μ2E424(α2r),A65=E423(β2r),A66=μ2E424(β2r)A67=−μ1E413(α3r),A68=−μ1E414(α3r),A69=−μ1E423(β3r),A610=−μ1E424(β3r)A71=A72=0,A73=μ2E723(α2r),A74=μ2E724(α2r),A75=E723(β2r),A76=μ2E724(β2r)A77=−μ1E713(α3r),A78=−μ1E714(α3r),A79=−μ1E723(β3r),A710=−μ1E724(β3r)A81=A82=0,A83=μ2E823(α2r),A84=μ2E824(α2r),A85=E123(β2r),A86=μ2E824(β2r)A87=−μ1E813(α3r),A88=−μ1E814(α3r),A89=−μ1E823(β3r),A810=−μ1E824(β3r)

{A91=A92=A93=A94=A95=A96=0,A97=E113(α3r)A98=E114(α3r),A99=E123(β3r),A910=E124(β3r)A101=A102=A103=A104=A105=A106=0,A107=E413(α3r)A108=E414(α3r),A109=E423(β3r),A1010=E424(β3r)

where: μ1=G2/G1,μ2=G3/G1, the superscript l = 1 denotes the Bessel function of the first kind, l = 3 represents the Hankel function of the first kind, and l = 4 indicates the Hankel function of the second kind.

References

1. Guan X, Yang N, Zhang W, Li M, Liu Z, Wang X, et al. Vibration response and failure modes analysis of the temporary support structure under blasting excavation of tunnels. Eng Fail Anal. 2022;136(15):106188. doi:10.1016/j.engfailanal.2022.106188. [Google Scholar] [CrossRef]

2. Li Y, Jin X, Lv Z, Dong J, Guo J. Deformation and mechanical characteristics of tunnel lining in tunnel intersection between subway station tunnel and construction tunnel. Tunn Undergr Space Technol. 2016;56:22–33. doi:10.1016/j.tust.2016.02.016. [Google Scholar] [CrossRef]

3. Pomasoncco-Najarro A, Trujillo-Valerio C, Arauzo-Gallardo L, Raymundo C, Quispe G, Dominguez F. Pre-split blasting design to reduce costs and improve safety in underground mining. Energy Rep. 2022;8(1):1208–25. doi:10.1016/j.egyr.2022.07.109. [Google Scholar] [CrossRef]

4. Lin Z, Jiang Y, Xu C, Chi M, Fang T, Guan L, et al. Performance of adjacent metro tunnels during deep excavation: a case study in Hangzhou. Transp Geotech. 2025;53:101585. doi:10.1016/j.trgeo.2025.101585. [Google Scholar] [CrossRef]

5. Fan Z, Zhang J, Xu H. Theoretical study of the dynamic response of a circular lined tunnel with an imperfect interface subjected to incident SV-waves. Comput Geotech. 2019;110:308–18. doi:10.1016/j.compgeo.2019.02.026. [Google Scholar] [CrossRef]

6. Pao YH, Mow CC, Achenbach JD. Diffraction of elastic waves and dynamic stress concentrations. J Appl Mech. 1973;40(4):213–9. doi:10.1115/1.3423178. [Google Scholar] [CrossRef]

7. Mow C, Mente L. Dynamic stresses and displacements around cylindrical discontinuities due to plane harmonic shear waves. J Appl Mech. 1963;30(4):598. doi:10.1115/1.3636625. [Google Scholar] [CrossRef]

8. Lu L, Yuan L, Tao G. Dynamic response sensitivity of urban tunnel structures under blasting seismic waves to parameters. Explos Shock Waves. 2014;34(6):701–8. doi:10.11883/1001-1455(2014)06-0701-08. [Google Scholar] [CrossRef]

9. Li XH, Long Y, Ji C, Zhou X, He YY, Lu L. Dynamic stress concentration factor for tunnel surrounding rock underblasting seismic waves. Rock Soil Mech. 2013;34(8):2218–24. doi:10.16285/j.rsm.2013.08.032. [Google Scholar] [CrossRef]

10. Lee VW, Trifunac MD. Response of tunnels to incident SH-waves. J Eng Mech Div. 1979;105(4):643–59. doi:10.1061/JMCEA3.000251. [Google Scholar] [CrossRef]

11. Xu H, Li T, Xu J, Wang Y. Dynamic response of underground circular lining tunnels subjected to incident P waves. Math Probl Eng. 2014;2014(1):1–11. doi:10.1155/2014/297424. [Google Scholar] [CrossRef]

12. Xu CJ, Ding HB, Tong LH, Luo WJ, Wang N. Scattering of a plane wave by shallow buried cylindrical lining in a poroelastic half-space. Appl Math Model. 2019;70:171–89. doi:10.1016/j.apm.2019.01.029. [Google Scholar] [CrossRef]

13. Tong L, Yu Y, Hu W, Shi Y, Xu C. On wave propagation characteristics in fluid saturated porous materials by a nonlocal Biot theory. J Sound Vib. 2016;379:106–18. doi:10.1016/j.jsv.2016.05.042. [Google Scholar] [CrossRef]

14. Ding H, Huang N, Muhammad C, Xu C, Tong L. Negative Poisson’s ratio locally resonant seismic metamaterials vibration isolation barrier. Acta Mech Sin. 2024;40(8):1–14. doi:10.1007/s10409-024-23370-x. [Google Scholar] [CrossRef]

15. Ding H, Huang N, Xu C, Xu Y, Cao Z, Zeng C, et al. A locally resonant metamaterial and its application in vibration isolation: experimental and numerical investigations. Earthq Eng Struct Dyn. 2024;53(13):4099–113. doi:10.1002/eqe.4214. [Google Scholar] [CrossRef]

16. Ding H, Yu Y, Xu C, Pu X, Guo W, Tong L. Analytical modeling for nonlinear seismic metasurfaces of saturated porous media. Int J Mech Sci. 2025;303:110666. doi:10.1016/j.ijmecsci.2025.110666. [Google Scholar] [CrossRef]

17. Xia Y, Han G, Mei W, Pan P-Z, Li M, Yan M, et al. Dynamic responses of a shallow lined tunnel with imperfect interface under transient P Wave. Int J Struct Stab Dyn. 2024;24(17):2450187. doi:10.1142/S0219455424501876. [Google Scholar] [CrossRef]

18. Tan Y, Yang M, Li X. Dynamic response of a circular lined tunnel with an imperfect interface embedded in the unsaturated poroelastic medium under P wave. Comput Geotech. 2020;122(1):103514. doi:10.1016/j.compgeo.2020.103514. [Google Scholar] [CrossRef]

19. Yi C, Lu W, Zhang P, Johansson D, Nyberg U. Effect of imperfect interface on the dynamic response of a circular lined tunnel impacted by plane P-waves. Tunn Undergr Space Technol. 2016;51(1):68–74. doi:10.1016/j.tust.2015.10.011. [Google Scholar] [CrossRef]

20. Fang XQ, Zhang TF, Li BL, Yuan RJ. Elastic-slip interface effect on dynamic stress around twin tunnels in soil medium subjected to blast waves. Comput Geotech. 2020;119(5):103301. doi:10.1016/j.compgeo.2019.103301. [Google Scholar] [CrossRef]

21. Xiang GL, Tao M, Zhao R, Zhao HT, Wu CQ, Memon MB. Dynamic characteristics of rockbolt anchorage structure under radial cylindrical P wave. Soil Dyn Earthq Eng. 2023;174:108176. doi:10.1016/j.soildyn.2023.108176. [Google Scholar] [CrossRef]

22. Peng S, Li X, Mitani Y, Gao J. Multiple-stage dynamic responses and failure behaviors of surrounding rocks subjected to development blasting: exponential and triangular paths. J Rock Mech Geotech Eng. 2025;17(6):3773–89. doi:10.1016/j.jrmge.2024.07.004. [Google Scholar] [CrossRef]

23. Cheng R, Chen W, Hao H, Li J. A state-of-the-art review of road tunnel subjected to blast loads. Tunn Undergr Space Technol. 2021;112(4):103911. doi:10.1016/j.tust.2021.103911. [Google Scholar] [CrossRef]

24. Li X, Li C, Cao W, Tao M. Dynamic stress concentration and energy evolution of deep-buried tunnels under blasting loads. Int J Rock Mech Min Sci. 2018;104(1):131–46. doi:10.1016/j.ijrmms.2018.02.018. [Google Scholar] [CrossRef]

25. Kavandi P, Ganjian N, Panji M. A DR-BEM approach for analyzing the transient SH-wave scattering problems: a comparative study. Eng Anal Bound Elem. 2024;169(3):105962. doi:10.1016/j.enganabound.2024.105962. [Google Scholar] [CrossRef]

26. Kavandi P, Panji M, Ganjian N, Marnani JA. An alternative dual reciprocity BEM for P-SV wave propagation problems: a comparative study. Eng Anal Bound Elem. 2025;176(A):106238. doi:10.1016/j.enganabound.2025.106238. [Google Scholar] [CrossRef]

27. Panji M. A half-space TD-BEM model for a seismic corrugated orthotropic stratum. Eng Anal Bound Elem. 2023;152:655–77. doi:10.1016/j.enganabound.2023.04.032. [Google Scholar] [CrossRef]

28. Mei W, Pan P, Wang Z, Xie Y, Feng Y, Chi F. Dynamic responses and failure characteristics of surrounding rock in a deep-buried tunnel with an arbitrary cross-section subjected to blasting P-wave. Comput Geotech. 2025;186:107447. doi:10.1016/j.compgeo.2025.107447. [Google Scholar] [CrossRef]