圆形洞室在径向非均匀荷载下的瞬态响应

耿大新， 陶 彪， 胡文韬

(华东交通大学 土木与建筑学院，南昌 330013)

随着地下空间深入开发，地下管线、地铁隧道、水下隧道等结构日益增多，弹性介质中含空腔或壳状结构的动力响应问题一直是研究的热点。然而这些结构多受一个随机的内部动态荷载作用，并且结构中大多数内源荷载并不是环向均匀的。有的结构由于受外部约束，荷载传递不均，向某方向集中汇聚。因此确定非均匀瞬态荷载引起的结构动力响应是地下工程领域中十分重要的问题。

针对圆形洞室与薄壁壳体等在轴对称瞬态荷载下的响应问题，目前主要的分析方法有解析法与数值法。解析法方面已有许多学者利用波函数展开法、积分变换法等方法进行研究。早在20世纪90年代初Senjuntichai等[1]采用波函数展开法，推导了全空间圆柱形腔体在三种不同类型轴对称荷载下的径向位移、应力、孔隙压力的精确通解，并通过数值反演拉普拉斯解，而得到时域解。在其基础上，Gao等[2-3]进一步研究了在衬砌洞室内部作用三种荷载情况下的动力响应。随后，Engin等[4]和陆建飞等[5]将洞室稳态下的解答，分别推广到半空间及任意洞室中，研究了弹性空间内土骨架的位移应力表达。此外，也有类似研究考虑了瞬态弹性波入射的情况，Karinski[6]、王滢[7]和李伟华[8]分别研究了衬砌洞室在不同种瞬态波散射下应力、位移、动应力集中的时域解，并考虑了刚度、衬砌厚度等因素对动应力集中的影响。翟朝娇等[9]针对反平面冲击荷载作用下洞室的瞬态响应问题进行探讨，研究分析了沿z轴方向瞬态荷载对土体动力响应的变化规律。然而当动荷载强度非常大，应变水平较高，对于这类问题一般常采用数值方法进行研究。Feldgun等[10]开创性的利用戈杜诺夫变分差分方法，分析了弹性、塑性及多孔介质中衬砌洞室的动力响应问题，并就解答的正确性与Glen等[11]的结果相印证。目前来说，已有研究大都针对洞室在环向均匀冲击荷载下的动力响应解析解求解[12-13]或者利用数值法结合有限元、边界元等对洞室内部均匀爆炸荷载作用下的动力响应分析[14]。但对径向非均匀均布荷载下的动力响应与波动特性，却一直未见类似研究，为此建立一种针对径向非均匀荷载作用下圆形洞室瞬态响应的计算方法，对隧道工程与地下空间领域具有深远意义。

本文将基于弹性介质波动理论，运用波函数展开法与Laplace变换法，根据圆形洞室内表面非均匀应力边界条件，求解出单位脉冲荷载下圆形洞室的动力响应解答。给出了全空间洞室中应力和位移场在时域内的数值解，并通过算例，分析了径向非均匀瞬态荷载下的波动特性以及剪切模量、不同角度对应力位移场的影响。

1 圆形洞室模型及波场求解

假定岩体为单相弹性介质，无限长圆柱形洞室埋置其中。因而，洞室内表面作用非对称瞬态荷载的动力响应问题可简化为平面应变问题。径向非均匀荷载随时间t,环向角度θ变化，如图1所示，r为极轴，a为圆形洞室内半径。

图1 径向非均匀荷载洞室模型Fig.1 Local concentrated load cavern model

图2 瞬态三角形脉冲荷载Fig.2 Transient triangular pulse load

2 岩体的控制方程

非均匀瞬态荷载F(t)从洞室内部传递至洞室边界后，在岩体中产生向外传播的膨胀波，洞室外的岩体视为单相弹性介质，其几何方程为[15]：

(1)

假定符合理想线弹性模型关系为：

σij=λ1δijεkk+2μ1εij

(2)

σr=(λ1+2μ1)εr+λ1εθ

(3a)

σθ=(λ1+2μ1)εθ+λ1εr

(3b)

σrθ=2μ1εrθ

(3c)

在极坐标系下，其振动方程可表示为：

(4)

(5)

式中：εr、εθ、εrθ分别为岩体径向应变、环向应变、切向应变；σr、σθ、σrθ为岩体径向应力、环向应力、切向应力；λ1为岩体的Lame常数；δij为Kronecker参数，当i≠j时δij=0，i=j时δji=1；ur、uθ为岩体介质的径向位移、环向位移。

将式(1)～(3)代入方程式(4)～(5)，得到以位移表示岩体的控制方程为：

(6)

(7)

由于ur(r,θ,t)和uθ(r,θ,t)是相互耦合的，为了解耦引入岩体部分的位移标量势函数φ(r,θ,t)和矢量势函数ψ(r,θ,t)，根据Helmholtz矢量分解定理有：

(8)

位移表示为：

(9)

对时间t进行Laplace变换和逆变换为：

(10)

(11a)

(11b)

由式(11)整理可得下式：

(12)

由式(12)可得到势函数Laplace变换后，满足如下的Helmholtz方程：

(13a)

(13b)

(14a)

(14b)

在极坐标下，Laplace算子与Laplace变换后的势函数可表示为：

(15a)

(15b)

采用分离变量法，对于线性系统中变换后势函数的解可表达为如下形式[16]：

(15c)

将式(15)代入式(13)，整理后可分解为两个不同的波数方程：

(16a)

(16b)

式(16)是n阶虚宗量Bessel函数，ki为波数，下标i=1,2。

式(16)中，势函数在极坐标下的通解可用Bessel函数线性组合的形式表达：

In(k1r)(A2cosnθ+B2sinnθ)]

(17a)

In(k2r)(C2cosnθ+D2sinnθ)]

(17b)

式中：Ai、Bi、Ci、Di，i=1,2为待定系数，In(·)为第一类虚宗量Bessel函数，Kn(·)为第二类虚宗量Bessel函数。

根据本文假设，在无限空间中，在r→∞时，须满足u,v→0，In不满足假设，因此待定系数A2、B2、C2、D2=0，则岩体势函数可表示为：

(18)

将式(18)代入表达式(9)，并考虑本构关系式(2)，可得极坐标下岩体中位移、应力以势函数的表达如下：

(19)

3 边界条件以及数值求解

本文研究在无限空间中半径为a的圆柱形洞室在内部受非均匀性冲击荷载作用如图1所示。确定待定系数与波场关系后，利用边界条件求解上述势函数中的待定系数，考虑圆形洞室与岩体交界面的边界条件可得：

(20)

本文为求径向非均匀瞬态荷载下的动力响应表达，荷载形式如下所示，将脉冲荷载表达式进行Laplace变换，得到Laplace变换域下表达，其表达式为：

(21)

式中：f(θ)可为任意形式的径向非均匀动荷载；b,c为径向非均匀荷载形状参数；T为三角形脉冲荷载的周期。

(22)

式中：P11、P12、P13、P14、E11、E12、E13、E14为系数项，具体表达详见附录。通过矩阵求解出待定系数后，代入式(19)中即可求出频域下应力与位移的解。

求得频域解答后，由于半解析解的形式直接进行Laplace逆变换较为困难，利用Laplace数值逆变换，转换为时域中的解，本文采取的是Durbin[17]数值逆变换方法，其变换表示为：

(23)

根据收敛准则确定NSUM的范围：

(24)

4 计算结果与算例分析

4.1 结果验算

(1)为了验证本文计算的合理性与正确性，将本文计算结果与文献[1]结果进行对比。为此将本文求解的径向非均匀瞬态荷载退化为径向均匀的瞬态荷载，即选取b=c=1，F0=0.1 MPa，λ=1.73×108Pa，a=3。计算所得无量纲环向应力、径向位移随时间分布图，如图3所示。本文在t*=8与t*=18峰值位置结果稍小，原因是参考文献[1]采用的数值逆变换方法与本文不同，在数值逆变换时无法完全拟合出对应的实部项α值。其次激励函数局部时间换算至全局时间过程取值偏大，将造成波数到达峰值时间略微推后。由图3可知，本文在环向应力的变化趋势上与文献[1]基本一致，说明了本文公式推导结果的合理性。

图3 本文退化为均匀的瞬态荷载与文献[1]结果比较Fig.3 Comparison of variation of hoop stress with time between present work and Ref.[1]

(2)另将本文计算结果与文献[18]数值模拟结果相对比，进一步验证推导结果的正确性。本文模型取值令F0=0.068 7 MPa，a=0.02 m，b=5,c=4，土体物理力学参数等与文献[18]一致，计算所得沿洞室径向不同位置处的正应力，如图4所示。图中r表示θ=0°位置上距洞室内表面的距离，由图可知正应力分布曲线与文献[18]数值模拟结果吻合较好，由此说明了本文推导结果的正确性。

图4 本文计算结果与文献[18]结果比较Fig.4 Result comparisons between calculation results presented in this paper and those in Ref.[18]

4.2 算例分析

考虑洞室r=a处(内表面)各个角度θ=0°,30°,60°,90°对土体应力位移响应的影响，土体的基本参数，见表1。

表1 计算参数Tab.1 Calculation parameters

图5表示洞室内壁，不同角度上的应力与位移响应值。从图5(a)曲线可知，t*>10时各角度位移响应曲线基本趋近于0。随着θ角0°→90°变化，环向位移先增大后减小，0°、90°的响应值为零。由于环向相互挤压变形，导致uθ随时间增大，各角度到达峰值的时间不同并且数值上差距较大。特别地，当接近t*=2时，30°处位移开始减小，60°位置处的位移刚到达最大值。表明在30°位置位移开始减小时，使得同一内径上60°位置的环向位移值有一小幅上升的阶段，每个θ角的振动是独立且异步的。图5(b)中可以看出不同角度对径向位移的影响十分显著。当θ从0°→30°时，径向位移峰值减小近20%。而θ持续增加到90°时，位移峰值衰减速率逐渐放缓，愈接近90°位置，衰减越慢，θ对径向位移的影响越小。与环向位移不同是各角度到达峰值的时间相同，振动是同步的。其中θ=0°时位移的幅值最大，与假设激励函数性质相同。在考虑径向非均匀的脉冲荷载最大位移响应时，应充分考虑0°即荷载集中位置的情况。

图5(c)、(d)应力在t*=1附近到达最大值，且各角度到达峰值的时间都相同，环向与径向均处于受压状态。随时间推移，环向应力由相互挤压状态变为环向拉伸状态，直至趋于稳定。θ角0°→90°时应力幅值逐渐减小，变化规律基本相同。并且θ改变时不仅影响应力峰值，也导致不同角度上应力衰减的速率有所差异。由图(c)、(d)可知θ从0～90°方向上其衰减速率逐渐减小，在θ=0°方向上能量扩散速度最快。图5(e)所示剪切应力随时间波动递减，呈往复态势。应力减小至t*=2时有一显著增大过程，越靠近90°位置往复性越明显，剪切响应越小，直到趋近0。此现象类似于环向应力，不同之处是环向应力随时间推移沿洞室环向拉压状态改变，而剪应力虽呈往复态势，但剪切方向不会改变。

图5 不同角度对洞室内表面(r=a)位移应力响应的影响Fig.5 Displacement stress response of the chamber surface at different angles

考虑不同径向距离对土体应力位移响应的影响，图6给出r*=r/a,r分别取1倍、1.2倍、1.4倍、和1.6倍洞径，角度取θ=0°或30°的应力位移随时间变化分布曲线。由于θ=0°时环向位移为零，无法考虑环向位移的波动关系，故任取一角度研究其波动关系，此算例取30°时的环向位移进行分析，下同。值得注意的是，从图(a)、(b)位移曲线来看，尽管r*增加位移响应进入峰值的时间各不相同，但是不同r*对应的曲线几乎是同一时间衰减至0。距离荷载中心位置越远，振动到达峰值的时间越长，振动衰减的越快。

从径向与环向应力图来看，图6(c)、(d)曲线都是随r*增大响应逐渐减小。径向应力波动周期比环向周期明显要短，在t*=2左右就衰减至零。环向应力波动时间之所以更长，是因为在同一土环内，环向会受到相互挤压作用，造成波动持续的时间更长。然而无限空间下，径向不存在相互作用，应力随着入射波不断向远处传播而越来越小，波动周期的也更短。

图6 不同径向距离对同一角度下位移应力响应的影响Fig.6 Displacement stress response at different radial distances at the same angle

对不同拉梅常数下的应力位移时域响应进行分析，剪切模量μ0=1.15×108Pa，μ取0.1、0.3、0.5倍的剪切模量，θ分别取0°或30°，其他参数如表1所示。图7(a)、(b)计算结果表明不同剪切模量对响应有显著影响。对于三角形脉冲荷载，剪切模量增加不仅使位移响应值提前到达最大值，并且也导致峰值大幅度降低。这表明介质刚度越大，对变形吸收越多，位移响应也逐渐减弱。而剪切模量变化对位移响应的衰减速率没有影响。当剪切模量较小时，位移响应的周期也越长，剪切模量越大，波动的时间越短，响应越早趋于平稳。

图7(c)为剪切模量对应力响应的影响，由图可知剪切模量的增加使环向应力值提前到达最大值，当剪切模量增大至0.3μ0后，剪切模量对响应到达峰值的时间的影响降低，0.3μ0以后的曲线几乎同一时刻到达峰值。而应力峰值的变化规律与位移响应规律相反，随剪切模量增大应力峰值越来越大，响应也愈加明显。从图中曲线可知，剪切模量对其衰减速率有显著影响，剪切模量增大时，应力衰减速率逐渐增大，相对于位移曲线没有表现出这种规律。并且位移与应力曲线都是在t*=1左右达到峰值，符合时域特征激励函数的假设。

图7 不同剪切模量对同一角度下位移应力响应的影响Fig.7 Displacement stress response of different shear modulus at the same angle

考虑不同时刻下的沿洞室环向应力位移时域响应，取洞室内表面r=a处。土体参数见表1，t*分别取0.5 s、1 s、3 s、5 s、10 s五个瞬时点。从图8(a)可知，在波动初始阶段t*=0.5时，环向在接近15°位置位移最大，在t*=1时环向位移最大值在30°。当t*=3时，环向位移最大值在60°位置出现，这说明随时间推移环向位移极值会随角度发生改变，各时刻下每个角度上的振动都是独立的，并不是同步到达峰值后衰减。在施加荷载初始阶段t*<3时，uθ出现峰值的角度随时间增大而增大。径向位移曲线分布如图(b)所示，不同t*下，最大位移出现在洞室θ=0°和θ=180°位置。时间变化并不影响径向位移峰值的位置，并且越接近90°位置位移响应衰减越多。从位移响应可知，在t*>3后位移波动趋于平稳，径向位移的响应值远大于环向位移。

图8(b)、(c)为环向应力与径向正应力沿洞室内表面的分布。在三角形脉冲荷载作用下，环向应力状态会发生改变，初始阶段环向受压，随环向响应增大，环向应力逐渐变为受拉。并且环向应力的极值随时间转动。当t*=0.5，环向应力在0°和180°取极大值，当t*=3时，极大值发生在90°和270°位置。与环向应力规律不同的是径向正应力绝大数都是受压状态，时间不会改变径向拉压状态。由此可以看出在洞室内表面位置，不同瞬时的差别，直接影响曲线出现极值的位置。

图8 不同瞬时对洞室内表面(r=a)位移应力响应的影响Fig.8 Displacement stress response of the chamber surface at different time

5 结 论

本文基于弹性动力学理论，建立无限介质中圆形隧道的径向非均匀瞬态荷载模型。采用波函数展开法，并利用三角函数正交性，求得洞室在非均匀瞬态荷载下频域内的半解析解，通过Laplace数值逆变换得到时域下应力位移的解答。此外，通过算例还分析了不同角度、介质模量等对响应的影响。并得到了以下结论：

(1)径向非均匀荷载作用下，环向位移响应随时间推移，各个角度的振动都是异步的，在各角度的振动到达峰值的时间各不相同，其中角度为0、π/2、π等位置处的环向位移为零。环向应力在拉压状态改变后峰值衰减也具有异步性，在径向应力和位移中并未出现类似现象。且在径向非均匀荷载作用时，由于环向响应的异步性，也使得结构发生破坏的几率增大。

(2)在径向非均匀荷载作用时，洞室环向应力与位移的极值位置随时间推移发生旋转。在t*<1荷载施加初始阶段，环向受到相互挤压，应力极值所在位置在0-π内随时间增大而增大。t*>1荷载释放过程中，环向由受压状态变为受拉状态，位移响应逐渐增强，位移的极值位置也开始旋转，位移旋转速率明显慢于应力。

(3)剪切模量对圆形洞室响应的影响较大，不仅使应力位移响应的峰值提前到达，并且极大的影响响应的幅值。剪切模量越大时，位移幅值越低，应力的幅值越高。

(4)环向位移与径向相比，其响应恢复到稳定的时间更长，环向振动的周期大于径向。由于应力拉压状态发生改变，当t*>4时，洞室内表面出现微小的负位移。θ在0与π荷载径向集中处，不论是径向或环向应力其波动幅值明显大于其他方向，且径向应力、位移的最值大于环向。

附录：