非圆形隧洞解析应力空间分布特征

2022-11-04杜小洲胡志平王振林安学旭陶磊

科学技术与工程 2022年27期

杜小洲，胡志平，王振林，安学旭*，陶磊

(1.陕西省引汉济渭工程建设有限公司，西安 710024； 2.长安大学建筑工程学院，西安 710061)

深埋高地应力条件下隧洞开挖经常伴随有围岩岩爆、剥落、大变形等一系列灾变现象，已给地下隧洞安全建设构成了重大挑战[1-2]，虽然引起这些灾变现象的因素较多，但现有理论和实践证明这些灾害与隧洞掌子面推进过程围岩重分布空间应力分布密切相关[3]，因此，有必要对深埋复杂应力场环境中隧洞围岩空间应力分布特征进行研究。

关于隧洞围岩应力分布特征，学者们采用了各种方法进行了大量研究。理论方面，已有诸多学者采用弹塑性理论对非均匀应力场中圆形隧洞断面应力分布特征及塑性区形状进行了各种研究[4-9]。对于复杂非圆形隧洞断面，Lü等[10]、李岩松等[11]、李元辉等[12]通过复变函数理论对复杂隧洞断面围岩应力分布特征进行了分析。目前虽然解析理论较多，但多集中在二维平面中。数值模拟方面，因其成本低、可重复以及较强的信息输出能力已被广泛地应用于去分析隧洞开挖过程中出现的各种工程问题，通常采用边界应力释放法、岩芯刚度折减法或岩芯置换技术等效手段来再现掌子面推进过程中的空间效应作用[13-14]。如陶志刚等[15]、蔡武强等[16]等通过三维数值模型分析了深埋岩体隧洞开挖过程中围岩应力、位移以及掌子面的非线性挤出效应。虽然数值模拟具有诸多优点，由于隧洞开挖是一个三维空间问题，对三维模型进行过细网格划分必然引起较高计算时间成本。试验方面，一些学者在制作物理试验模型时，通过在预定位置预埋与开挖洞型一致的柱体，等试件干燥后将其拔出或顶出，然后再进行加载，或在预定位置通过手动或利用小型机械设备开挖出一定形状的隧洞来模拟并分析掌子面推进过程中围岩应力空间分布特征[17-18]。三维模型试验研究隧洞掌子面空间效应虽效果较好，但其缺陷也较显著，如造价高、可重复利用性差等。

由上述分析可知，对于隧洞围岩应力分布特征的研究虽然方法较多，但各有局限。现基于收敛约束法，将非圆形隧洞围岩应力及位移解析理论与可视化软件Tecplot通过程序相结合，提出一种高效、直观的隧洞解析应力空间分布方法，并采用该方法对引汉济渭工程岭北隧洞开挖过程中掌子面附近围岩空间应力分布特征以及不同初始主应力场水平倾角下的隧洞围岩空间应力分布特征进行研究。

1 隧洞解析应力空间分布方法

隧洞开挖过程中由于开挖面的卸荷作用使得围岩中产生应力重分布现象并伴随位移释放，且整个过程不会在开挖瞬间立即完成，而是在掌子面推进过程中逐渐完成，如果在无支护条件下，当掌子面推进一定距离后整个位移释放才完成[19]。目前，这个过程可通过收敛约束法中纵向变形剖面(longitudinal deformation profile,LDP)与围岩特征(surrounding rock characteristic，GRC)曲线来表示[20]，纵向变形曲线能反映与掌子面不同距离处隧洞断面围岩位移，一般可通过实测数据或数值模拟等方法获取。围岩特征曲线反映了隧洞断面围岩位移与虚拟支护力之间的关系[21]，本文中采用隧洞断面位移释放系数与虚拟支护力的之间的关系来反映，该曲线可通过考虑位移释放系数的隧洞围岩应力及位移解析理论获得。

解析应力空间分布模型构建：首先采用相同的变量位移释放系数，建立隧洞LDP曲线与GRC曲线之间的关系如图1所示，然后将平面空间位移释放系数对应至三维空间中支护与掌子面之间的距离，确定三维隧洞空间任意点处围岩应力及位移值，最后通过可视化软件Tecplot输出三维空间应力分布图，可供实际工程支护设计与施工提供参考。

为了使该方法通用化，具有可复制性，将解析理论通过编程在MATLAB软件中实现快速计算，并将围岩力学计算结果数据以及围岩纵向变形曲线数据导入可视化软件Tecplot，即可输出三维隧洞空间应力图，整个计算过程只需要几分钟即可完成，计算流程图如图2所示。

η(x)=ur/urmax, ur为掌子面附近任意隧洞断面围岩位移，urmax为无支护条件下隧洞断面最大围岩位移；Pi为虚拟支护力，MPa; x(i)为距离掌子面的距离，m图1 隧洞LDP与GRC曲线之间的关系Fig.1 Relationship between LDP and GRC curves

α为初始主应力场水平倾角，(°)；σ1和σ3分别为最大和最小主应力，MPa；E1和E2分别为围岩与衬砌的杨氏模量，MPa；υ1和υ2分别为围岩与衬砌的泊松比图2 隧洞围岩解析应力空间分布方法流程图Fig.2 Flow chart of analytical stress spatial distribution method for tunnel surrounding rock

2 考虑位移释放系数的隧洞解析方程

2.1 隧洞围岩位移解析表达

对于如图3(a)所示的隧洞断面竖向主轴与其主应力场非平行的隧洞，可通过三角函数分解为如图3(b)所示的力学模型。

σy为初始主应力场水平应力分量，MPa；σx为初始主应力场竖向应力分量；τxy为初始主应力场剪切分量，MPa图3 非圆形隧洞围岩应力场Fig.3 Surrounding rock stress field of a non-circular tunnel

L21表示z平面中支护外边界，L11表示z平面中支护内边界；L22表示映射ζ平面中支护外边界，L12表示映射ζ平面中支护内边界；θ为映射ζ平面中极坐标角度,(°)；ρ为映射ζ平面中极坐标半径，mm；Ro为映射ζ平面中支护内边界半径，mm；i为复数图4 z平面上隧洞断面形状映射为ζ平面上单位圆周Fig.4 Tunnel section shape on the z plane is mapped to the unit circumference on the ζ plane

如图4所示，对于单连通域内复杂隧洞断面形状可通过保角映射转换为单位圆外域，映射函数可以表示为

(1)

式(1)中：ω(ζ)为关于变量ζ的函数；R为反映隧洞大小的系数；ck为反映隧洞形状的系数。

实际工程中，隧洞断面虽然形态较多，但大多关于竖轴对称，因此，式(1)中映射函数系数ck可用常实数表示。

对于无支护隧洞，围岩应力求解需通过两个解析函数获得，即

φ1(ζ)=Γ1ω(ζ)+φo(ζ)

(2)

ψ1(ζ)=(Γ2+iτxy)ω(ζ)+ψo(ζ)

(3)

式中：Г1=σx(1+λ)/4，Г2=σx(λ-1)/2；λ=σx/σy；φo(ζ)和ψo(ζ)分别为圆外解析函数，φo(ζ)、ψo(ζ)。

可以表示为

(4)

(5)

(6)

式中：ɑk(k=1,2,3,…)为解析函数的系数；j=2,3,…,n-2；系数lk求解见参考文献[10]。

由于模型应力边界条件中存在非对称初始主应力场剪应力分量，因此式(4)解析函数中系数ɑk一定是复数。ɑk可通过线性方程组求解，即

Xikak=Bi，i,k=1,2,…,n-2

(7)

式(7)中：当i=k时，Xik=kLi+k+1，若i+k+1>n，Xik=0；当i≠k时，Xik=kL2i+1-1，若i+k+1>n，Xik=-1。Bi列矩阵为

(8)

式(8)中：i=2,3,…,n-2。

由式(1)、式(4)及式(5)可知无支护隧洞围岩中任意点处位移为

(9)

式(9)中：μ1和ν1分别为竖向和水平总位移，mm；G1=E1/[2(1+υ1)]；κ1=3-4υ1；E1和υ1分别为围岩杨氏模量与泊松比。

隧洞开挖过程中随着掌子面的推进，围岩中会产生应力重分布，围岩中位移会逐渐释放，在支护前假定隧洞围岩产生的位移为

μ2+iv2=η(μ1+iv1)

(10)

式(10)中：η为位移释放系数，0≤η≤1.0；μ2和ν2分别为支护前围岩中产生的竖向与水平位移，mm。

支护施工后，由于围岩位移的进一步释放使围岩与支护之间会产生相互作用，由于支护约束作用在围岩中产生的位移可表示为

(11)

式(11)中：μ3和ν3分别为支护约束作用在围岩中产生的竖向和水平位移，mm；Φ2(ζ)和ψ2(ζ)分别为支护对围岩的约束作用的圆外解析函数，可表示为

(12)

(13)

式中：bo、bk、do、dk均为待求复系数。考虑到支护对隧洞无穷远处的影响较小，因此，认为在无穷远支护对围岩的作用为0，由此可得k1bo=k1(bo1+ibo2)=conj(do)= (do1-ido2)，conj(i)表示对系数i共轭。

支护在围岩作用下任意点处位移可表示为

(14)

式(14)中：G2=E2/[2(1+υ2)],κ2=3-4υ2，E2和υ2分别为支护结构的杨氏模量和泊松比；μ4和ν4分别为支护中产生的竖向和水平位移，mm；Φ3(ζ)和ψ3(ζ)分别为支护在围岩作用下的圆外解析函数，可表示为

(15)

(16)

式中：po、ek、fk、qo、gk、hk均为待求复系数。

2.2 边界连续条件

为了使得解析函数满足应力及位移边界连续条件，则解析函数[式(12)、式 (13)、式(15)、式(16)]在支护外边界L21处应符合边界应力及位移连续条件，在支护内边界L11处应符合边界应力连续条件。

在支护外边界L21处由隧洞径向应力连续条件可得

(17)

式(17)中：σ=eiθ。

在支护边界L21处由隧洞径向位移连续条件可知

(μ1+iv1)(1-η)+(μ3+iv3)=μ4+iv4

(18)

将式(9)、式(10)、式(11)、式(14)代入式(18)化解可得

(19)

将ζ平面中支护内边界用Roσ表示，由支护内边界处L11处应力连续条件可知

(20)

将式(17)、式(19)及式(20)左右两侧ζ±j(j=0,1,…,N)系数相同的实部和虚部整理联立方程可获得12N+6个方程和12N+6个未知量，这些线性方程组可统一表示为

Aijxj=Fi

(21)

式(21)中:Aij为方程组左侧项系数矩阵；Fi为方程组右侧项列矩阵，整个方程组在MATLAB中通过程序进行求解。

2.3围岩应力及位移求解

2.3.1 支护

将式(21)求解的变量代入式(12)、式(13)、式(15)、式(16)解析函数，结合式(1)可得

σpL+σθL=Re[φ′|3(ζ)/ω′(ζ)]

(22)

(23)

(24)

式中：σθL、σρL、τρθL分别为支护环向、径向及剪切应力，MPa。

2.3.2 围岩

同理将式(21)求解的变量代入式(13)和式(14)，可得围岩中应力计算公式为

σpR+σθR=4Re[φ′|4(ζ)/ω′(ζ)]

(25)

(26)

φ4(ζ)=φ1(ζ)+φ2(ζ)

(27)

ψ4(ζ)=ψ1(ζ)+ψ2(ζ)

(28)

(29)

式中：σθR、σρR、τρθR分别为围岩中环向、径向及剪切应力，MPa。

上述解析理论中，N越大所求得的计算结果越精确，实际上N>50基本可以满足多数工程计算精度要求，本文中N=120。

3 算列分析

3.1 工程概况

引汉济渭工程是陕西境内重大水利设施项目，隧洞全长81.779 m，其中秦岭岭北硬岩段隧洞地质环境比较突出，埋深最深达2 012 m，岩石平均抗压强度为170 MPa，隧洞开挖过程中发生岩爆、大变形、岩体剥落、卡机等诸多工程灾变问题。为了给隧洞设计及施工提供依据，在岭北6号实验洞附近进行了地应力测试，测试结果如表1所示[22]。围岩杨氏模量与泊松比参数分别采用E1=64.5 GPa，υ1=0.19，支护杨氏模量与泊松比参数分别采用E2=30 GPa，υ2=0.20。隧洞断面形式如图5所示，断面映射函数通过复合形优化算法获得，具体公式为

z=ω(ζ)=3.983 77(ζ-0.051 9+

0.001 7ζ-1+0.039 1ζ-2-0.042 2ζ-3+

0.022 9ζ-4-0.001 1ζ-5-0.007 1ζ-6)

(30)

表1 6#试验洞地应力测试结果Table 1 Situ stress test results of No. 6 test tunnel

图5 引汉济渭工程隧洞断面Fig.5 Tunnel section of Hanjiang to Weihe River valley water diversion project

3.2 隧洞断面空间应力分布特征

根据3.1节隧洞断面形状及基本力学参数，通过数值模拟法获得了隧洞顶部位移释放系数与掌子面不同距离之间的关系曲线如图6所示。

将该曲线与解析理论相结合，应用第1节提出的隧洞解析应力空间分布法可获得引汉济渭工程隧洞断面空间应力分布云图如图7所示。

如图7所示，隧洞掌子面推进过程中，隧洞左拱脚及右拱肩处环向应力集中较大，尤其是左拱脚，最大可达5.0倍的初始主应力场竖向应力分量，而右拱脚及左拱肩处环向应力集中相对较小。沿隧洞轴线纵向，距离掌子面0.5D范围内，隧洞左拱脚处环向应力梯度变化较大，右拱肩处变化其次，而左拱肩处变化相对较小。距离掌子面0.5～2.0D范围内，隧洞左拱脚及右拱肩处环向应力集中区域随掌子面的推进在逐渐增大。距离掌子面2.0D范围外，隧洞环向应力分布区域基本趋于稳定。由此可见，沿隧洞轴纵向，逐渐远离掌子面时，隧洞环向集中应力分布呈现出3阶段特征，应力集中程度快速增长区、应力集中范围扩展区、应力集中稳定区。

图6 隧洞顶部位移释放系数与掌子面之间的关系Fig.6 Relationship between displacement release coefficient of tunnel top and tunnel working face

D为隧洞断面最大洞径图7 隧洞空间应力分布云图Fig.7 Spatial stress distribution of tunnel

4 初始主应力场水平倾角的影响

为了探究隧洞初始主应力场水平倾角α对隧洞解析应力空间分布特征的影响，设置了如表2所示的3种应力边界工况。

基于3.1节隧洞断面形状及力学参数，通过数值模拟获得了不同主应力场水平倾角条件下的隧洞顶部位移释放系数与掌子面不同距离之间关系曲线如图8所示，相应曲线方程如表3所示。

将表3中拟合曲线方程分别与隧洞解析理论相结合，通过应用第1节提出的隧洞解析应力空间分布方法可获得不同初始主应力场条件下隧洞围岩解析应力空间分布云图如图9所示。

表2 模型边界不同应力场参数Table 2 The different stress field parameters of the model boundary

图8 隧洞顶部位移释放系数与掌子面之间的关系Fig.8 Relationship between displacement release coefficient of tunnel top and tunnel working face

表3 拟合曲线方程Table 3 Fitting curve equations

图9 不同初始主应力场水平倾角下隧洞环向应力集中空间分布云图Fig.9 Spatial distribution nephogram of circumferential stress concentration of tunnel under different horizontal inclination angles of initial principal stress field

如图9所示，在α=0° 时，隧洞空间应力集中区域主要分布在隧洞左右拱脚及拱肩处；在α=30° 时，隧洞空间应力集中区域分布主要集中隧洞左拱脚及右拱肩处，而隧洞左拱肩及右拱脚处环向应力集中程度与α=0° 时相比，在降低；在α=60° 时，隧洞空间应力集中区域分布主要表现在左拱脚及拱顶偏右侧，与α=30° 时相比，隧洞顶部偏右侧环向应力集中程度在升高。由此可见，初始主应力场水平倾角变化对隧洞环向应力环向空间分布位置影响较大。沿隧洞轴线纵向逐渐远离掌子面时，不同初始主应力场水平倾角下隧洞环向应力集中区域也呈现出应力集中梯度快速升高、应力集中区域逐渐扩展、应力集中区域稳定三阶段特征，说明初始主应力场水平倾角变化对隧洞环向应力集中纵向变化影响较小。