浅埋地下孔洞、夹杂或隧道等地下结构对地震波散射在抗震和抗爆研究中具有重要的理论意义和应用价值，也一直是土木工程和地震工程领域中备受关注的热点问题^[1-5]。为揭示浅埋地下结构对地震波散射的机理，学者们采用了各种数学方法开展了大量的解析求解工作。比如，基于波函数展开法、Graf公式和镜像方法，Lee等^[1-2]最早推导了弹性半空间中圆形孔洞和隧道对平面SH波散射的闭合解，并分析了衬砌模量、隧道埋深、入射波频率和角度等因素对隧道周边动力响应的影响。采用类似方法，袁晓铭^[5]推导了地下圆形夹杂对平面SH波的散射的闭合解，分析了夹杂模量和埋深等因素对地面运动的影响。刘殿魁等^[6]和王国庆等^[7]采用复变函数法给出了地下单个和多个圆形孔洞对平面SH波的散射问题的级数解，计算了孔洞周边的动应力集中系数。陈志刚等^[8-9]结合保角变换和复变函数法求解了各向同性和各向异性弹性半空间内任意形状孔洞对平面SH波的散射问题。李敏等^[10]和刘刚等^[11]采用复变函数方法分别求解了半圆形凸起地形与地下圆形孔洞的组合以及等腰三角形凸起地形与地下圆形孔洞的组合对平面SH波的散射问题。齐辉等^[12]基于复变函数方法推导了半圆形凸起地形与地下圆形孔洞的组合对柱面SH波散射的级数解。Gao等^[13]借助分区方法和多极坐标转换技术推导了地下马蹄形孔洞对平面和柱面SH波散射问题的级数解。Chen等^[14]基于波函数展开法推导V形河谷地形与地下圆形孔洞、夹杂和隧道的组合对平面SH波散射问题的级数解。Zhang等^[15]基于复变函数方法推导了各向异性半空间内非完全粘结隧道对平面SH波散射问题的级数解。对于入射SH波作用下的反平面运动问题，只涉及单个标量波动方程，求解相对简单。而对于入射P波、SV波或Rayleigh面波等作用下的面内运动问题，由于耦合了纵波波速和横波波速，难以通过分离变量的方法直接构造出满足地震波半空间表面边界条件的波函数，给问题的求解带来很大的困难。通过将半空间表面假定为一个很大的圆弧（即“大圆弧假定”），Lee等^[16]推导出弹性半空间内圆形孔洞对平面SV波散射的近似级数解。随后，梁建文等^[17-18]基于“大圆弧假定”推导了弹性半空间内洞室群对平面P波和SV波散射问题的近似解析解。Mei等^[19]结合波函数展开法和积分变换，推导了瞬态平面P波作用下地下圆形孔洞动力响应问题的近似解析解。Lin等^[20]用Hankel函数的波数积分形式构造出适应地面边界条件的波函数，得到了弹性半空间内圆形隧道对平面P波散射的精确解。Liu等^[21]利用保角变换将地表和孔洞表面映射成一个圆环，推导出弹性半空间内圆形隧道对平面P波、SV波和Rayleigh面波散射的解析解，与边界元解^[22]对比具有较好的一致性。随后，Liu等^[23]进一步给出了饱和多孔弹性半空间内圆形隧道对平面P1波和SV波散射的解析解。

尽管已经有许多关于浅埋地下结构对地震波散射问题的解析解，但大多数是针对平面波且衬砌与围岩完全粘结的情况。只有少数考虑了波源距离^[12-13]和衬砌与围岩非完全粘结^[15]的影响。通常情况下，天然地震波的波源一般距离地下结构较远，波阵面的半径远大于地下结构，波线到达地下结构时接近互相平行，因此可将入射波考虑为平面波。然而，人工地震波（如爆破产生的）的波源通常距离地下结构较近，波阵面与地下结构尺寸差异不大，波线到达地下结构时的角度差异较大，因此不能将入射波简单视为平面波。一般的方法认为衬砌和围岩是完全粘结接触，但由于微裂纹和间隙介质，衬砌和围岩并不总是完全粘结接触的^[15]，而围岩与衬砌的接触条件与衬砌内的位移和应力密切相关^[24]。为研究波源距离和非完全接触的综合效应对隧道动力响应的影响，并出于简便考虑，本文采用反平面线源表示入射波，采用位移不连续模型描述衬砌和围岩的不完全粘结界面；基于波函数展开法、Graf公式和镜像方法，推导反平面线源荷载作用下浅埋非完全粘结圆形隧道动力响应问题的级数解，旨在揭示波源距离和粘结程度等因素对隧道动力响应的影响规律，以期为地下结构抗震和抗爆设计提供一定参考依据。

1.   理论模型

反平面线源荷载作用下浅埋圆形非完全粘结隧道动力响应的理论模型如图1所示。围岩和衬砌均为各向同性线弹性介质，围岩的剪切模量、密度、剪切波速分别为μ₁、ρ₁和c₁，衬砌的剪切模量、密度和剪切波速分别为μ₂、ρ₂和c₂。衬砌的外径和内径分别为a₁和a₂，隧道中心的埋深为H，线源的埋深为D，线源与隧道中心的水平距离为L。采用位移不连续模型表示衬砌的与围岩的接触条件，用接触刚度K_s表示粘结程度。为便于描述波场，在隧道中心、线源以及隧道中心正上方地表分别建立笛卡尔坐标系(x_1, y₁)、(x_2, y₂)和(x_, y)，以及极坐标系(r_1, θ₁)、(r_2, θ₂)和(r_, θ)。线源在极坐标(r_, θ)中的位置可表示为(r_0, θ₀)。

图  1  理论模型示意图

Figure  1.  Schematic diagram of theoretical model

下载: 全尺寸图片 幻灯片

上述问题的控制方程为^[25]：

式中：下标 j =1,2，分别表示围岩和衬砌；w_j表示z方向的位移；β_j= ω/c_j为波数，ω为圆频率；δ(·)表示狄拉克δ函数。

地表的应力自由边界条件为：

式中：${\tau }_{\theta {\textit{z}}}^{\left(1\right)}$为围岩中的周向剪应力。

衬砌内表面的应力自由边界条件为：

式中：$\tau_{r_1{\textit{z}}_{1}}^{(2)}$为衬砌中的径向剪应力。

衬砌外表面与围岩的接触条件为：

式中：w⁽¹⁾和w ⁽²⁾分别为围岩和衬砌中的位移，${\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(1\right)}$为围岩中的径向剪应力。

如图1所示，围岩中包括线源激发的入射波w^(inc)，地表激发的反射波w^(ref)，以及衬砌外表面激发的外行散射波w^(a)和相应地表激发的外行散射波w^(b)。衬砌中只包含驻波w^(c)。

首先，线源激发的入射波w^(inc)的数学形式可表示为^[25]：

式中：i为虚数单位，${\mathrm{H}}_{0}^{\left(1\right)}\left(\cdot \right)$表示第零阶第一种汉克尔函数。

为便于表示反射波w^(ref)，在地表上方建立线源的虚像，并在虚线源上设置笛卡尔坐标系(x_4, y₄)和极坐标系(r_4, θ₄)，可得w^(ref)的数学形式：

定义自由场位移w^(ff)=w^(inc)+w^(ref)，并利用坐标关系式：x₄=2D−x₂，y₄=y₂，x₂=x+D,y₂=y+L，x=x₁−H，y=y₁，可得w^(ff)在极坐标系(r₁, θ₁)下的数学表达式：

式中：${R}_{\mathrm{1,2}}={r}_{1}^{2}\pm 2\left(D\mp H\right){r}_{1}\mathrm{cos}{\theta }_{1}+2L{r}_{1}\mathrm{sin}{\theta }_{1}+{\left(D\mp H\right)}^{2}+{L}^{2}$。

相应的应力${\tau }_{{r}_{1}{{\textit{z}}}_{1}}^{\left(\mathrm{f}\mathrm{f}\right)}$为：

式中：${\mathrm{H}}_{0}'^{\left(1\right)}(\cdot )$表示${\mathrm{H}}_{0}^{\left(1\right)}\left(\cdot \right)$的一阶导数。

将式(7)和(8)在[−π, π]上分别对角度θ₁进行傅里叶级数展开，可得：

其中：

衬砌外表面激发的外行散射波w^(a)的数学表达式可写为：

式中：A_n和B_n为待定系数。

相应的应力${{\tau }}_{r_1{\textit{z}}_{1}}^{\left(\mathrm{a}\right)}$和$\tau _{\theta_1{\textit{z}}_{1}}^{\left(\mathrm{a}\right)}$为：

根据镜像方法，可假定地表上存在一个虚拟隧道，如图1所示，其激发的外行散射波w^(b)可表示为：

相应的应力${{\tau }}_{{{r}}_{3}{\textit{z}}_{3}}^{\left(\mathrm{b}\right)}$和${{\tau }}_{{{\theta }}_{3}{\textit{z}}_{3}}^{\left(\mathrm{b}\right)}$为：

由于w^(b)数学形式基于极坐标系(r₃, θ₃)，而边界条件式(3)和(4)基于极坐标(r₁, θ₁)，因此应将w^(b)变换到极坐标(r₁, θ₁)上，一种可行的方法是运用Graf公式^[5]：

式中：J_m(·)是第m阶贝塞尔函数；当m=0时，ε_m=1，当m≠0时，ε_m=2。

结合式(16)和(19)可得w^(b)在极坐标(r₁, θ₁)下数学形式：

式中：${U}_{m,n}={\varepsilon }_{m}\left[{\mathrm{H}}_{m+n}^{\left(1\right)}\left(2{\beta }_{1}H\right)+{\left(-1\right)}^{n}{\mathrm{H}}_{m-n}^{\left(1\right)}\left(2{\beta }_{1}H\right)\right]/2$，${V}_{m,n}={\varepsilon }_{m}\left[{\mathrm{H}}_{m+n}^{\left(1\right)}\left(2{\beta }_{1}H\right)-{\left(-1\right)}^{n}{\mathrm{H}}_{m-n}^{\left(1\right)}\left(2{\beta }_{1}H\right)\right]/2$。

相应的应力${{\tau }}_{{{r}}_{1}{\textit{z}}_{1}}^{\left(\mathrm{b}\right)}$为：

衬砌包含内外两个表面，因此衬砌中的散射波w^(c)为驻波场，其数学形式可表示为：

式中：C_n, D_n, E_n和F_n为待定系数。

相应的应力${\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(\mathrm{c}\right)}$和${\mathrm{\tau }}_{{\mathrm{\theta }}_{1}{{{\textit{z}}}}_{1}}^{\left(\mathrm{c}\right)}$为：

至此，已完全构建围岩和衬砌中的波场，其中围岩中的总位移波场w⁽¹⁾=w^(ff)+w^(a)+w^(b)，衬砌中的总位移波场w⁽²⁾=w^(c) 。

下面将利用衬砌内外表面的边界条件式(3)和(4)建立求解待定系数A_n、B_n、C_n、D_n、E_n、F_n的代数方程组。

首先将等式${\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(2\right)}{=\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(\mathrm{c}\right)}$和式(23)代入边界条件式(3)，并依据三角函数的正交性，可得如下方程组：

接着将等式w⁽¹⁾=w^(ff)+w^(a)+w^(b)，w⁽²⁾=w^(c)，${\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(1\right)}{=\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(\mathrm{f}\mathrm{f}\right)}+{\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(\mathrm{s}\right)}$和${\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(2\right)}={\mathrm{\tau }}_{{{r}}_{1}{{{\textit{z}}}}_{1}}^{\left(\mathrm{c}\right)}$和式(9)～(10)、(13)～(14)和(22)～(23)代入边界条件式(4)，并依据三角函数的正交性，可得如下方程组：

联立式(25)、(27)和(29)可得如下方程组：

式中：m = 0, 1, 2···。

同理，联立式(26)、(28)和(30)可得如下方程组：

式中：m = 1,2,3···

通过截断和求解方程组(31)～(33)和方程组(36)～(38)，可以计算出待定系数A_m、B_m、C_m、D_m、E_m、F_m。接着将这些系数代入式(13)、(16)和(22)中，就可以确定围岩和衬砌中的散射场。再结合自由场式(7)，可以求出整个模型内的波场。

为了缩减变量数目，定义如下无量纲频率：

式中：f为入射波频率，λ₁=c₁/f为入射波波长。

2.   精度校验

为校验本文级数解的精度，定义衬砌内边界应力残差以及衬砌外边界位移和应力残差如下：

式中：${w}_{0}={w}^{\left(\mathrm{inc}\right)}\left({\beta }_{1}\sqrt{{\left(D-H\right)}^{2}+{L}^{2}}\right)$，τ₀=μ₁β₁W₀。

不失一般性，本节模型参数取值为：ρ₁= ρ₂ =2500 kg/m³，μ₁ =10.4 GPa，μ₂/μ₁ = 1.5，a₁=5 m，a₂/a₁=0.9，H/a₁=2.0，D/a₁=5.0和L/a₁=5.0。图2和图3分别给出了当η=1.0、K_s/μ₁=1.0 m⁻¹和η = 3.0、K_s/μ₁ = 0.1 m⁻¹时，${\Delta }{{\tau }}_{{{r}}_{1}{\textit{z}}_{1}}^{\left({{a}}_{2}\right)}$、$\mathrm{\Delta }{{w}}^{\left({{a}}_{1}\right)}$和${\Delta }{{\tau }}_{{{r}}_{1}{\textit{z}}_{1}}^{\left({{a}}_{1}\right)}$随级数解的截断项数N的关系。由图2(a)和图3(a)可知，在截断项数较小的情况下，衬砌内边界的应力残余量极小，几乎可以忽略，但随着截断项数的增加应力残余量有增大趋势。再由图2(b)和2(c)以及图3(b)和3(c)可知，在截断项数较小的情况下，衬砌外边界的位移余量和应力余量相对较大，不过随截断项数的增加，这些数值会显著减小。另外，比较图2和图3的结果可发现，随着η的增加，达到同等精度水平级数解的截断项数也要增加。比如，当η=1.0时，N=7可保证足够精度（所有边界条件残余量小于1%），而当η=3.0时，需要N=13才能保证同等精度。综合上述结果可以推论，本文所提出的级数解是正确且收敛的，但应适当选取截断项数N，以保证解的精度。

图  2  当η = 1.0和K_s/μ₁ = 1.0 m⁻¹时衬砌内外边界的位移和应力残差与截断项数N之间的关系

Figure  2.  The relationship between the displacement and stress residuals of the inner and outer boundaries of the lining and the number of truncated terms N when η=1.0 and K_s/μ₁=1.0 m⁻¹

下载: 全尺寸图片 幻灯片

图  3  当η = 3.0和K_s/μ₁ = 0.1 m⁻¹时衬砌内外边界的位移和应力残差与截断项数N之间的关系

Figure  3.  The relationship between the displacement and stress residuals of the inner and outer boundaries of the lining and the number of truncated terms N when η=3.0 and K_s/μ₁=0.1 m⁻¹

下载: 全尺寸图片 幻灯片

3.   结果和讨论

通过对上述级数解进行参数分析，本节将系统讨论衬砌和围岩的接触刚度、衬砌模量、衬砌厚度、隧道埋深和线源距离等因素对隧道动力响应的影响规律。为了精简内容，本节部分参数固定取值为：η=1.0、a₁=5 m、μ₁=10.4 GPa、ρ₁=ρ₂=2500 kg/m³。

3.1   衬砌和围岩的接触刚度对隧道动力响应的影响

为分析衬砌和围岩的接触刚度对隧道动力响应的影响，图4给出了当μ₂/μ₁=1.5、a₂/a₁=0.9、D=H=2a₁、L/a₁=5.0时，不同接触刚度K_s/μ₁情况下衬砌内表面的归一化位移和周向剪应力幅值的分布。可以看出，衬砌和围岩的接触刚度对隧道的动力响应具有较大影响。当接触刚度非常小时（K_s/μ₁=0.01 m⁻¹，对应黑色实线），衬砌内表面的位移和周向剪应力幅值相对入射波幅值基本可以忽略。说明此时入射波的大部分能量被衬砌散射回围岩，只有极少部分能量进入衬砌。当接触刚度较大时（K_s/μ₁=1.0或10 m⁻¹，分别对应蓝色和黑色的虚线），衬砌内表面的位移和周向剪应力均较大，且K_s/μ₁=1.0 m⁻¹和K_s/μ₁=10 m⁻¹的结果差异很小。这说明当接触刚度达到K_s/μ₁=1.0 m⁻¹时，衬砌与围岩已接近完全粘结。当接触刚度较小时（K_s/μ₁=0.1 m⁻¹，对应红色的实线），衬砌内表面的位移和周向剪应力呈现出异常高的幅值，显著大于接触刚度较大（K_s/μ₁=1.0或10 m⁻¹）的情况，且背波侧（0°<θ₁<180°）的幅值并未低于迎波侧（180°<θ₁<360°）。这可能是因为地震波在衬砌内出现了多重散射，散射波的多次叠加导致衬砌出现整体振动而产生的现象。以上结果表明，要提高隧道减震效果，必须控制好衬砌与围岩的接触刚度。

图  4  衬砌内表面归一化位移和周向剪应力与衬砌和围岩的接触刚度的关系

Figure  4.  The relationship between normalized displacement and circumferential shear stress on the inner surface of lining and the contact stiffness between lining and surrounding rock

下载: 全尺寸图片 幻灯片

3.2   衬砌模量对隧道动力响应的影响

除了衬砌与围岩的接触刚度外，衬砌模量的大小对隧道的动力响应也有较大影响。图5给出了当K_s/μ₁=1.0 m⁻¹、a₂/a₁=0.9、D=H=2a₁、L/a₁=5.0时，不同衬砌模量μ_2/μ₁情况下衬砌内表面的归一化位移和周向剪应力幅值的分布。由图可知，随着衬砌模量增大，衬砌内表面的位移幅值会减小，但周向剪应力幅值会增大，且周向剪应力的增大程度要大于位移的减小程度。因此，只增大衬砌模量并不能有效地提高隧道的抗震性能。

图  5  衬砌内表面归一化位移和周向剪应力与衬砌模量的关系

Figure  5.  The relationship between normalized displacement and circumferential shear stress on the inner surface of lining and lining modulus

下载: 全尺寸图片 幻灯片

3.3   衬砌厚度对隧道动力响应的影响

相较于衬砌模量，衬砌厚度更容易进行调整。图6给出了当K_s/μ₁=1.0 m⁻¹、μ_2/μ₁=1.5、D=H=2a₁、L/a₁=5.0时，不同衬砌内外径比a₂/a₁（该数值越大，衬砌厚度越小）情况下的衬砌内表面归一化位移和周向剪应力幅值的分布。由图可知，在其余条件相同情况下，增大衬砌厚度能同时降低衬砌内表面的位移和周向剪应力。因此，增大衬砌厚度可有效提高隧道的抗震性能。但考虑到隧道的开挖成本和空间使用要求，适度增大衬砌厚度并优化减震层的设计（即控制围岩和衬砌的接触刚度），可能是更经济的手段。

图  6  衬砌内表面归一化位移和周向剪应力与衬砌厚度的关系

Figure  6.  The relationship between normalized displacement and circumferential shear stress on the inner surface of lining and lining thickness

下载: 全尺寸图片 幻灯片

3.4   隧道埋深对隧道动力响应的影响

隧道埋深对地震波的传播路径具有显著的影响，从而影响隧道的动力响应。图7给出了当K_s/μ₁= 1.0 m⁻¹，a₂/a₁ = 0.9，μ_2/μ₁ = 1.5，D/a₁ = 2.0 、L/a₁ = 5.0时，不同隧道埋深H/a₁情况下的衬砌内表面归一化位移和周向剪应力幅值的分布。首先，由图7(a)可知，当隧道埋深很浅时（H/a₁=1.5，对应黑色实线），衬砌内表面的最大位移出现在衬砌的左下侧。随着隧道埋深增加，最大位移幅值出现的位置会逐渐向上移动，当隧道埋深达到H/a₁=3.0时（对应黑色的虚线），最大位移出现隧道拱顶的左侧附近。这是因为当道埋深H/a₁=1.5时，线源在隧道中心左下方，由于入射波直达衬砌左下侧，因此此处位移较大。但随着隧道埋深增加，线源会逐渐移动到隧道中心左上方，因此衬砌左上方位移会逐渐增大。接着，从图7(b)可看出，随着隧道埋深增加，衬砌内表面的最大周向剪应力也会逐渐上移，但一直保持在隧道左下侧范围内。此外，随隧道埋深增加周向剪应力总体上呈现减小趋势。综上可知，在其他因素不变的情况下，增大隧道埋深有助于提高隧道的抗震性能。

图  7  衬砌内表面归一化位移和周向剪应力与隧道埋深的关系

Figure  7.  The relationship between normalized displacement and circumferential shear stress on the inner surface of lining and tunnel burial depth

下载: 全尺寸图片 幻灯片

3.5   线源距离对隧道动力响应的影响

线源与隧道的距离会影响入射波波阵面的曲率和射线方向，从而对隧道动力响应产生影响。图8给出了当K_s/μ₁=1.0 m⁻¹、a₂/a₁=0.9、μ_2/μ₁=1.5、D=H=2.0a₁时，不同线源水平距离L/a₁情况下衬砌内表面的归一化位移和周向剪应力幅值的分布。从图中可看出，线源距离对衬砌内表面的位移和周向剪应力的分布有很大的影响。在线源距离较小的情况下（L/a₁= 2.0，对应黑色实线），衬砌的迎波侧（180°<θ₁<360°）位移和周向剪应力幅值显著大于背波侧（0°<θ₁<180°）。这是因为当线源距离隧道较近时，入射波的大部分能量都被衬砌朝反方向散射，只有少部分能量通过衍射进入背波侧。随着线源距离增加，入射波的曲率减小，射线方向逐渐趋近水平，衬砌对入射波的屏障作用减弱，因此衬砌的背波侧位移和周向剪应力幅值逐渐增大。另外，由于地面反射增强的原因，又以拱顶附近的增加程度最为显著。上述结果说明，对隧道进行抗震设计时，还应考虑波源距离的影响。

图  8  衬砌内表面归一化位移和周向剪应力与线源距离的关系

Figure  8.  The relationship between normalized displacement and circumferential shear stress on the inner surface of the lining and the distance from the line source

下载: 全尺寸图片 幻灯片

4.   结　论

采用位移不连续模型表示衬砌和围岩的接触条件，首先建立了反平面线源荷载作用下浅埋圆形非完全粘结隧道动力响应的理论模型，然后通过波函数展开法对改模型进行了求解了，接着通过衬砌内外表面的边界条件校验了级数解的精度和收敛性，最后通过参数分析系统讨论了衬砌与围岩的接触刚度、衬砌模量、衬砌厚度、隧道埋深和线源距离等因素对隧道动力响应的影响，并得出以下结论：

(1)当接触刚度非常小时（K_s/μ₁=0.01 m⁻¹），衬砌和围岩接触面具有很好的隔振作用，隧道的动力响应较小；当接触刚度较大时（K_s/μ₁=1.0 m⁻¹），衬砌与围岩的接触界面接近完全粘结；当接触刚度为较小时（K_s/μ₁=0.1 m⁻¹），衬砌内表面的位移和周向剪应力呈现异常高的幅值，且衬砌左右两侧的幅值差别较小；

(2)增大衬砌模量会减小衬砌内表面的位移幅值，但同时会增大衬砌内表面的周向剪应力幅值；

(3)增大衬砌厚度能同时减小衬砌内表面的位移和周向剪应力幅值；

(4)在线源位置不变的情况下，增大隧道埋深会使衬砌内表面的最大位移和周向剪应力幅值向上移动，且周向剪应力幅值整体上趋于减小；

(5)在隧道埋深不变的情况下，增大线源与隧道的水平距离会使衬砌拱顶和背波侧的位移和周向剪应力相对幅值增大。