考虑矿柱流变效应的采空区力学模型分析

蔡 泳 邓国梁 甄利兵 成 墙 黄金涌

（1.煤矿灾害动力学与控制国家重点实验室,重庆 400044;2.重庆大学资源与安全学院,重庆 400044;3.贵州紫金矿业股份有限公司,贵州 贞丰 562200）

采矿方法中,房柱法划分矿块为矿柱和矿房,回采矿房而预留下矿柱,形成大量的采空区[1]。如果不及时处理采空区,可能造成人员上的伤亡和财产损失[2-5]。因此采空区自稳时间的预估显得尤其重要。在稳定时间内,对采空区进行处理,如支护、充填等工作,可避免安全事故的发生。采空区力学分析能有效地预估采空区稳定时间。多数力学分析将矿柱简化为winkler弹性基础和蠕变模型[6-8],将顶板简化成弹性薄板、弹性岩梁或弹性厚板组建力学模型分析[9-11]。Zhao等[12]基于岩梁模型,确定了顶板岩块转角与挠度的关系。孙琦等[13]将采空区顶板简化成弹性薄板,将矿柱简化成Kelvin-Voigt模型,考虑水平应力和垂直应力的共同作用,建立起多向荷载作用下的采空区力学模型,经过理论计算与现场实际监测结果对比分析,验证了模型的实用性。楼小明等[14]将矿柱简化成Hooke-Kelvin-kelvin蠕变模型,顶板简化为弹性薄板,建立力学模型应用到工程中,并指出矿柱的流变变形是采空区失稳的主要原因。Wang等[15]基于Reissner厚板理论建立采空区力学模型,通过互等定理推导出在不同的边界条件顶板挠度随时间变化的表达式,借助数值模拟软件验证解析解的合理性。谢学斌等[16]将采空区顶板和矿柱简化成弹性薄板和Burgers蠕变模型,构建出三维蠕变模型,借助尖点突变理论对采空区失稳倾向做出定量分析,并分析了各个影响因子对系统稳定性的影响。上述矿柱—顶板力学模型建立过程中,大多是将矿柱视为流变模型,顶板视为弹性模型。实际上,采空区自稳时间是受矿柱和顶板共同流变作用的影响。若对采空区影响因素进行敏感性分析,对主控因素考虑其流变性作用,对非主控因素理想化考虑,会缩减力学模型计算量的同时,还能保证模型精确度。

鉴于此,本研究以某钨矿2线0沿东采空区为工程背景,根据因素敏感性分析结果,构建采空区三维流变力学模型,得出采空区自稳时间的表达式,并利用Flac3D对采空区蠕变过程中应力、位移进行了分析。

1 采空区影响因素敏感性分析

1.1 工程背景

某钨矿采用房柱法开采。由于多年的开采和早年开采不规范,已经形成了大大小小的采空区,并且出现了矿柱片落、劈裂和顶板冒顶的现象。迫切需要对某钨矿2线0沿东采空区进行长期稳定性分析,采空区示意图如图1所示。

图1 2线0沿东采空区示意Fig.1 Schematic diagram of goaf along E0 of line 2

取某钨矿采场基本参数,分析采场的长×宽为80 m×80 m,矿柱高度6 m,埋深200 m。矿柱设计为5 m×5 m方柱,矿房跨度10 m。

1.2 影响因素敏感性分析

参考李令鑫等[17]对采空区安全系数的定义,将采空区的影响确定为6因素5水平,进行正交极差分析,所得结果如表1所示。

表1 各因素极差分析结果Table 1 Range analysis results of various factors

通过极差分析可得各个影响因素的主次关系:矿柱的宽高比＞矿柱的单轴抗压强度＞顶板上覆岩层的容重＞采空区埋深＞矿房跨度＞顶板厚度。

2 采空区长期稳定性力学分析

采空区的稳定与其结构参数和时间有关,考虑时间效应的作用下,矿柱能承受最大应力值应为岩石的长期强度[18]。根据因素敏感性分析结果,在采空区自稳时间的力学模型构建中,矿柱若简单视为弹性元件,可能会对结果造成较大的偏差。同时,顶板厚度影响性最小,力学模型的构建中可以将顶板简化为弹性薄板。

2.1 矿柱蠕变模型

经过文献[19]研究,对比不同应力下的蠕变实验曲线与蠕变模型发现Burgers模型能够较好表征硬质矿柱的蠕变行为。Burgers流变模型如图2所示。

图2 岩石Burgers流变模型Fig.2 Burgers creep model of rock

Burgers蠕变模型本构方程:

式中,η1、η2为牛顿黏性系数;k1、k2为弹性元件的弹性系数;σ为模型的总应力,MPa;ε为模型的总应变。

化简有:

2.2 采空区力学模型

根据因素敏感性分析,采空区的力学模型可以简化矿柱为Burgers模型,而顶板可以简化成弹性薄板。简化后的采空区力学模型如图3所示。

图3 采空区简化模型Fig.3 simplified model of goaf

随着开采的进行,顶板边界逐渐破坏,采空区的边界条件分为3个阶段:固支边、简支边和自由边。

（1）阶段一:固支边位移方程。当顶板处于固支边时,则顶板的控制方程为:

式中,D为顶板抗弯刚度,kN/m2;w为顶板挠度,m;ξ为将矿柱中的集中力均布化系数;σ为矿柱中应力,MPa;q为上覆岩层作用在顶板上的压力,MPa;n为采场内矿柱个数;A为矿柱的横截面积,m2;a、b为矩形采场的长边和短边,m;E为顶板岩石弹性模量,GPa;h为顶板厚度,m。

联立式（1）、式（3）消除式中σ,得到采空区顶板—矿柱系统的控制方程。

对式（4）,采用伽辽金法求其弱解。阶段一时。顶板边界处于固定,挠度和转角为0:

构造近似解的形式,设为

式中,w0（t）为顶板中心点随时间变化的挠度,m;φ（x,y）为挠度随位置变化的分布函数。当t=0时,固支边条件下采场的挠度如图4所示。

图4 t=0固支边时采场挠度Fig.4 Stope deflection with fixed support when t=0

近似解满足边界条件,加权余数表示为

式中,w0为采场中心挠度,m;J1、J2、J3为化简后系数。

（2）阶段二:简支边位移方程。阶段一采场位移达到破坏条件,此时,顶板边界出现了转角,处于简支边,边界的挠度和弯矩为0。

阶段二的挠度方程形式设为

简支边条件下的分布函数φ（x,y）如图5所示。

图5 t=0简支边时采场挠度Fig.5 Stope deflection with simply supported edge when t=0

同理,阶段二采场中心挠度的微分方程为

顶板继续持续破坏,表面裂隙形成并贯通,顶板已经完全失去了承载能力（D=0）。顶板的边界处于自由边,阶段三状态。

2.3 采空区自稳时间解析解推导

采空区的失稳是一个逐渐发展的过程,自稳时间可以看作阶段一、阶段二的时间总和。

采空区处在阶段一时,采场的挠度表示为

表达式满足边界条件,代入某钨矿基本数据和矿柱流变数据（表2）并根据式（7）有:

表2 矿柱流变参数Table 2 Rheological parameters of pillars

问题转化成解此微分方程的解,式（13）解得:

接收现场数据，并对数据进行分析，根据特定的计算方式给出相应的指令，这一部分称之为运算部分。逻辑控制部分：与现场设备相连接，接收运算部分的动作指令，判断现场设备的实际运行状况并将新的动作指令分别下发给各终端，使其按照要求单独动作或者联动动作。这2个部分互有区别，不应混为一谈。

式中,A1、A2为待定系数,根据边界条件确定。

t=0时,根据前人的研究[20]和Burgers蠕变方程,采场矿柱刚受压时,有初始位移和初始蠕变速度。

式中,σ为矿柱中应力,MPa。

求解出待定系数,最终得出阶段一时采场中心挠度随时间的表达式:

阶段一的破坏条件为

解得阶段一采空区自稳时间为32.14个月。

同理,也可以得到阶段二顶板中心挠度的表达式和破坏表达式。

式中,[σt]为允许抗拉强度,MPa,ν为泊松比。

代入数据,阶段二的破坏时间为88.03个月。则采空区总的自稳时间120.22个月。根据阶段一和阶段二的表达式,顶板中心点的挠度随时间变化的趋势如图6所示。

图6 顶板中心挠度随时间变化的曲线Fig.6 Curve of central deflection of roof with time

3 采空区长期稳定性数值分析

Flac3D是ltasca公司研发的连续介质力学分析软件[21],被大量用于岩土工程和采矿工程中。建立计算模型如图7所示,矿柱设为Burgers本构模型,顶板和底板为Mohr-Coulomb模型。设立5个矿柱顶部监测点,观察采空区失稳的破坏规律。蠕变计算一步设为150 h,计算至失稳后停止计算。

图7 Flac3D计算模型示意Fig.7 Schematic diagram of Flac3D model calculation

3.1 应力分析

开挖矿体对地下初始应力场造成扰动,导致地下的应力场重新分布。由图8和图9知,模型竖直应力云图上,竖直应力整体上呈现随着深度的增加而增加。采空区形成后,矿体采空部分形成了应力释放,矿柱中心部分有明显的应力集中。矿柱顶端的压应力为11 MPa接近理论计算的应力值。矿柱上部形成明显的拱形等值线,矿柱顶板中间出现了较小的拉应力。而随着蠕变迭代步数的进行,观察到矿柱上方的竖直拉应力和影响区域都在不断增加。

图8 开采完竖直应力云图Fig.8 Vertical stress nephogram after mining

图9 蠕变过程中竖直应力云图Fig.9 Vertical stress nephogram in creep process

由图10知,模型最大主应力云图上,采空区顶板上主要分布拉应力,最大拉应力出现在2个矿柱之间的中心的顶板上,未达到岩石的抗拉强度,此时顶板关键层并不会发生破坏。单独通过拉应力值的大小判断顶板是否会发生破断并不合理,还需要结合顶板下沉量综合分析。

由图11知,模型最小主应力云图上,最小主应力基本上呈层状分布,矿柱支撑位置出现了应力集中,矿柱支撑力分布表现出中间矿柱支撑力大,而边缘矿柱支撑力小的现象。

图11 蠕变过程中最小主应力云图Fig.11 Minimum principal stress nephogram in creep process

3.2 位移分析

由图12知,模型竖直位移云图上,矿房开采引起的应力变化相互叠加,导致周围岩体产生的位移也随之变化。此时,采空区顶板最大下沉量为18.93 mm,小于临界破坏条件21.30 mm。而采空区底板出现了正向位移,最大位移量为2.59 mm。在采空区的上部形成了规律的位移等值线,并随着埋深的增加而等值线的大小逐渐减小。矿柱的位移规律为:矿柱上半部表现为向下的位移,而下半部为向上的位移,表现出两端受压的状态。

图12 蠕变过程竖直位移云图Fig.12 Vertical displacement nephogram of creep process

对5个矿柱顶部中心位置监测位移,由图13知,1#监测点和5#监测点、2#监测点和4#监测点符合对称的结果。从整体上来看,采场中心的挠度到采场边缘呈现递减的现象,符合之前假设的势函数形式。5个监测点在开采完27.70个月内,矿岩在上覆岩层持续应力作用下显现出显著的时间效应,先变形速率快,矿柱的顶部位移一直呈上升趋势,后变形速率减小,曲线趋于平稳。最后在118.06个月时,模型停止计算。

图13 5个监测点竖直位移随时间的曲线Fig.13 Curve of vertical displacement of five monitoring points with time

4 结 论

（1）采用敏感性分析方法得出了各个影响因素的敏感性,根据敏感性分析结果将矿柱简化Burgers蠕变模型,顶板简化为弹性薄板,建立了采空区三维蠕变力学模型,结合Burgers蠕变方程,推导出顶板位移随时间变化的表达式。

（2）根据新构建的力学模型对某钨矿采空区自稳时间进行了预测,与数值模拟软件对比分析,验证了模型的合理性和可靠性。

（3）利用Flac3D软件对某钨矿采场模拟分析,随着蠕变迭代步数的进行,观察到矿柱上方的竖直拉应力和影响区域都在不断增加,顶板的位移变化速率呈先增大后减小的趋势。