基于复变函数理论的巷道围岩稳定性分析

2022-06-06王志强苏泽华

煤炭科学技术 2022年4期

王志强，黄鑫，武超，石磊，苏泽华

(1.中国矿业大学(北京) 能源与矿业学院，北京 100083；2.中国矿业大学(北京) 共伴生能源精准开采北京市重点实验室，北京 100083；3.中国矿业大学(北京) 煤炭安全开采与地质保障国家级实验教学示范中心，北京 100083)

0 引言

我国煤矿已进入深部开采阶段，部分矿井开采深度已超千米，深矿井开采，因其地压大，巷道经受较大的变形和破坏，造成的人员伤亡和经济损失也随之增加[1-2]，围岩塑性区的确定是巷道支护参数选取以及巷道稳定性评估的重要依据。孙广忠[3]结合弹性解构造应力分量表达式，在轴对称应力场塑性区表达式的基础上，得到围岩塑性区的范围；赵志强等[4-7]提出巷道“蝶形”破坏理论,为巷道围岩塑性区的分布规律以及巷道支护提供了有力的依据；文献[8-11]研究了两向不等压作用下圆形巷道围岩的弹塑性问题，在库仑屈服条件下,用摄动法给出弹塑性交界线的解析表达式；吕爱钟等[12]基于Mohr-Coulomb强度准则,解决了圆形隧洞在两向不等压地应力作用下的塑性区确定问题；张承客等[13]对硐室围岩弹塑性边界两侧分别采用复变函数理论和滑移线场理论得出弹性区和塑性区应力组合表达式,根据弹塑性边界上应力相等得到边界方程；袁文伯等[14]根据岩体应变软化的变形特性,建立了弹塑性软化模型，分析了软岩巷道围岩塑性区和破碎区的力学性态，推导出比卡斯特纳(Kastner)公式更符合实际、适用性更广泛的塑性区半径计算公式。文献[15-21]利用复变函数方法获得了Mohr-Coulomb强度准则下非轴对称弹塑性问题的近似解；董海龙等[22]以侧压系数为1的特例为基准,将依据近似隐式法得到的围岩塑性区半径与基于成熟轴对称平面应变理论的相应解析解进行对比；巷道围岩控制是矿井发展的基础，大部分学者对双向不等压圆形巷道围岩塑性区进行研究，但出于成巷难易程度考虑，很多矿井采用直墙半圆拱形巷道，鲜有学者对其塑性边界解问题进行相应的研究与探讨。

基于此，笔者以数值模拟结果为参考，采用复变函数理论，对直墙半圆拱形巷道围岩塑性区范围展开研究，以期为类似矿井巷道的理论研究和实践提供帮助。

1 巷道围岩应力复变函数解

复变函数解法在平面问题中的孔口问题中最有优越性，其特点在于通过保角变换，根据映射函数将复杂孔口(例如巷道断面常采用的拱形等非圆曲线)在z平面上所占区域，变换成ζ平面上所谓“单位圆”，即孔口边界变换为单位圆周界，有利于简化边界条件，从而推导出孔口围岩应力的弹性解析式。由于ζ平面为极坐标系，因此需将应力的复变函数表达式转化为ζ的极坐标函数。

根据弹性力学的极坐标应力变换公式得：

(1)

其中，η为z平面ρ轴与x轴的夹角，由式(1)得出：

(2)

令z=ω(ζ)，ζ为ζ平面内单位圆边界上任意一点，即ρ=1，ζ=ρ(cosθ+isinθ)=ρeiθ=eiθ，其中ρ、θ为点ζ关于坐标原点ζ=0的极坐标，则：

cos 2η+isin 2η=e2iη=

(3)

即得到ζ平面极坐标下的应力分量的复变函数表达式：

(4)

(5)

由弹性理论及复变函数理论可得：

(6)

(7)

令σ=ω(ζ)=eiθ，即σ为复变量ζ在巷道边界的值，则得到φ0(ζ)、ψ0(ζ)表达的边界条件，经柯西积分公式计算后得到：

(8)

(9)

f0为σ的已知函数，其表达式为：

(10)

建立双向不等压直墙半圆拱形巷道受力模型如图1所示。假设水平与竖直原岩应力均匀分布，巷道z轴方向长度远大于其他两方向长度，受力状态不予考虑，且巷道的埋深足够深，将此问题简化为无限区域内的平面应变问题。

图1 直墙半圆拱形巷道力学模型

其保角映射函数为

(11)

其中对于本文问题，R，ci都为实常数，正是要确定的量，在求解具体问题时，ci只能为有限值，在本文研究中i的最大值取4。

结合本文研究具体问题，联立式(6)、式(7)、式(8)，应用Harnack定理得：

(12)

将式(12)代入式(5)，可得：

(13)

cj1=-2δjcoskπ，

c1=c12+c22+c32-c11c21-(c11+c21)c31，

c2=c13+c23+c33+c11c22+

(c12-c11c21+c22)c31+c12c21+(c11+c21)c32，

c3=c14+c24+c34-c11c23-(c13+c23-c11c22+c12c21)c31+c12c22-c13c21+

(c21-c11c21+c22)c32-(c11+c21)c33

c4=c15+c25+c35+(c14+c24-c11c23+

c12c21-c13c21)c31+c12c22+c13c22+c14c21+

(c13+c23+c11c22+c12c21)c32+c11c24+

(c21-c11c21+c22)c33+(c11+c21)c34

令映射后的单位圆ρ=1，则ζ=eiθ=cosθ+isinθ，对于式(2)而言，由于巷道壁上σρ=0，则得σθ=4Re[Φ(ζ)]，代入式(13)，可以得出：

(14)

令

则得：

(15)

其中，

(16)

综合式(11)、式(15)、式(16)，最终得到直墙半圆拱巷道周边的切向应力计算公式为

σθ=4ReΦ(ζ)=

(17)

2 巷道围岩稳定性力学分析

2.1 巷道围岩极限平衡区计算

巷道开挖之后，围岩应力重新分布，巷道内部出现应力集中现象。其中当应力集中小于围岩极限强度时，巷道围岩处于弹性状态；而当应力集中大于极限强度时，围岩处于塑形状态，甚至处于破裂松动状态。巷道处于双向不等压前提下，铅垂方向巷道承载为σy，水平方向巷道承载为λσy，巷道承载如图2所示。

图2 巷道应力状态

根据极限平衡区内的应力平衡方程，可得：

(18)

其中，

(19)

式中，C为围岩黏聚力；φ为围岩内摩擦角。

将式(19)代入(18)得到：

(20)

对上式进行积分，得到：

(21)

(22)

当r=R时，σr=0。代入上式，得：

(23)

将式(23)代回式(22)，得到：

(24)

(25)

结合式(17)和(25)可得：

(26)

通过上述分析可知，直墙拱形巷道围岩的极限平衡区受到巷道宽度(2b)、直墙高度(H)、拱高(2h-H)、垂直应力(σy)、侧压系数(λ)、黏聚力(C)以及内摩擦角(φ)的影响。其中巷道围岩的黏聚力和内摩擦角可以通过巷道支护的方式来改变，也鉴于篇幅有限，此文仅对黏聚力和内摩擦角的大小对巷道围岩塑性区范围的影响进行计算研究，为类似工程提供理论支持。

2.2 巷道围岩塑性区计算

2.2.1 黏聚力对塑性区的影响

计算参数取为：σy=10 MPa，a=1.5 m，λ=1.8，k=0.3，h=2 m，H=3.2 m。内摩擦角取为30°，黏聚力分别取为1，3，5，7 MPa。运用Mathcad对式(26)代入上述参数进行计算，采用单因素分析法改变黏聚力的大小所获得的塑性区大小和形状如图3所示。

图3 不同黏聚力塑性区大小和形状

由图3可以看出，黏聚力越小，塑性区范围越大，并且在不同地应力方向，巷道围岩塑性破坏范围的增幅有明显的不同，在σy方向上的增大幅度较大，通过计算还可以发现，巷道顶板的塑性区范围略小于巷道底板，而巷道两帮的塑性区范围呈对称分布的趋势。

2.2.2 内摩擦角对塑性区的影响

黏聚力取为1.5 MPa，内摩擦角分别取为20°，30°，40°，50°，其他参数同2.2.1节)。运用Mathcad对式(26)代入参数进行计算，所获得的塑性区大小和形状如图4所示。

图4 不同内摩擦角塑性区大小和形状

由图4可以看出，内摩擦角越大，塑性区的范围越小，当内摩擦角比较小时，塑性区呈“X”形破坏，内摩擦角比较大时，塑性区呈类似于椭圆形破坏，随着内摩擦角的增大，巷道围岩塑性破坏范围的减幅有所减小，与黏聚力类似的是，塑性区范围在σy方向上的增大幅度较大。

3 算法的准确性分析

3.1 数值模拟分析

为验证算法的准确性，首先采用有限元软件ANSYS建立三维模型，使用FLAC3D数值模拟软件对直墙半圆拱形巷道周围塑性区分布进行数值模拟计算，整个模型的长×宽×高(X，Y，Z)为100 m×10 m×68.4 m，如图5所示。

图5 数值模拟模型

整体模型由223 440个单元和248 039个节点组成，网格大小从巷道向外逐步增加，模拟采用摩尔-库伦准则进行计算，根据表1所示的力学参数对各分组进行赋参。施加相应竖直向下的压力σy加载于模型顶部模拟未建覆岩重量，模型底部约束横向和纵向位移，前后左右约束横向位移。黏聚力和内摩擦角按照2.2节进行赋参，采用单因素分析法对模型分别赋以不同的黏聚力和内摩擦角的大小，研究其对巷道塑性破坏范围的影响。