动载扰动下深井护巷煤柱失稳破坏机理

邢闯闯 ，王 俊 ，宁建国 ，沈 圳

（1.山东科技大学 矿山灾害预防控制省部共建国家重点实验室培育基地, 山东 青岛 266590；2.山东科技大学 能源与矿业工程学院, 山东 青岛 266590）

0 引 言

近年来，煤柱型冲击地压灾害事故频率不断增加，严重威胁矿井的生产安全[1-2]。例如2016 年8月15 日0 时33 分，山东能源肥矿集团梁宝寺能源有限责任公司发生了煤柱型冲击地压事故，造成2人死亡。2019 年8 月2 日12 时24 分，开滦(集团)有限责任公司唐山矿业分公司也发生了煤柱型冲击地压事故，造成了12 人伤亡，直接经济损失高达614.024 万元。实践表明，煤柱型冲击地压一般具有隐蔽性强、无明显前兆特征等特点，造成设备破坏、巷道变形和人员伤亡等事故，给矿井的安全生产带来严重威胁。因此，有必要开展煤柱动力失稳机理研究，以期对煤柱型冲击地压预警与防治提供理论依据。

国内学者针对煤柱稳定性问题进行了大量的研究。部分学者从静力学角度进行了大量的试验。殷志强等[3]通过分离式霍普金森压杆试验，分析了不同静载下煤样的动态力学特性，得出峰值应力随静载增大而增大，峰值应变随静载增大而减小的结论。何永琛[4]通过对煤岩组合体试件进行单轴压缩试验，获得了不同煤岩组合体试件的破坏形态及峰值应变等变化特征。艾迪昊等[5]通过对不同粒度的型煤进行单轴压缩试验，得出型煤的破坏可分为压实，初始破裂，加速破裂和卸压4 个阶段。李成武等[6]通过室内试验和现场试验研究了煤体静载破坏中低频磁场时、频谱特征，并结合微震信号提出了低频磁信号的产生机制。

也有学者从动力学的角度，探讨了煤柱稳定性问题。徐青云等[7]采用理论分析、数值模拟和现场实测等方法，指出窄煤柱在受到剧烈采动影响时，主要在其中上部区域发生变形破坏。姜学伟[8]通过分析房采区煤柱在采动影响下围岩受力、煤柱破坏变形状态等得出了煤柱保持稳定性所必须的弹性区宽度计算公式，并作为判断煤柱是否破坏失稳的依据。李峰等[9]基于LS-DYNA 有限元分析软件的HJC 本构模型，模拟不同动载强度下煤柱的动态损伤过程。结果表明，应变率相同时，煤柱受动载荷越大，破坏越严重。王方田等[10]将动载荷以应力波的形式作用在煤体上，来分析不同动载荷烈度下煤柱坝体的动力响应机制，发现巷道围岩塑性区发育面积及深度显著增加，竖向应力积聚强度增大，并与动载荷烈度成正比。丛利等[11]通过数值模拟对相邻工作面回采中和回采后3 次高强度开采扰动对煤柱内应力变化情况的研究，得出煤柱内布置的联络巷会加剧巷道围岩的高应力集中，在外部动载扰动下，极易达到冲击载荷并诱发冲击的结论。因煤柱在外部动载扰动下的破坏失稳机理极其复杂，当前的理论研究和实践仍需进一步完善。因此针对不同尺寸深井护巷煤柱在外部动载扰动下的破坏特征和突变失稳机制开展了详细研究。

首先，采用数值模拟分析了动载冲击下煤柱失稳特征；其次，建立了动载作用下高静载煤柱突变失稳力学模型，推导了煤柱稳定性判别准则，揭示了动载作用下深井煤柱突变失稳机理，以期对煤柱稳定性控制提供参考。

1 深井护巷煤柱数值模型的建立

1.1 模型建立与参数选择

利用FLAC3D软件建立“顶板-煤柱-底板”数值计算模型，模型宽高比分别设置为4、6、8 和10，煤层高度保持为3 m。图1 为宽高比为10 的数值计算模型，模型共计9 层，7 种岩性，模型尺寸为34 m×34 m×39 m，共划分62 820 个单元网格。参考相关数值计算经验[12-13]，顶底板岩层采用各向同性弹性模型，煤层采用应变软化模型。煤层节理、裂隙发育，室内试验获得的煤样物理力学参数很难直接用于数值计算。因此，采用试错法拟合煤柱强度公式以校核煤柱力学参数。相关学者提出了多种煤柱强度公式[14]，其中Salamon 公式能够准确描述宽高比为2～20 时的煤柱强度，计算结果精确度能够满足现场要求，可用于计算现场保护煤柱的强度[15]。因此，采用Salamon 公式校正数值模型计算结果。

图1 数值计算模型Fig.1 Numerical calculation model

式中：Sp为保护煤柱强度，MPa；W为保护煤柱宽度，m；H为保护煤柱高度，m。

采用试错法时，首先建立“顶板-煤柱-底板”数值计算模型；其次，在模型顶部施加速度边界条件以模拟单轴压缩，获得煤柱单轴压缩强度；再者，通过调整煤层各项力学参数，使数值计算获得煤柱单轴压缩强度与Salamon 公式匹配。图2 为数值计算获得的煤柱强度与Salamon 公式对比图，由图2 可知，数值计算获得的煤柱强度与Salamon 公式基本保持一致，因此，采用表1 所示的岩性参数模拟煤柱的应变软化行为。

图2 数值模拟和Salamon 公式计算峰值强度对比Fig.2 Comparison of peak strength between numerical simulation and Salamon formula calculation

表1 煤柱模拟参数Table 1 Simulation parameters of protective coal pillar

1.2 边界条件及模型方案

以新河煤矿6302 工作面南部轨道集中巷，6301运输巷和6302 轨道巷之间的保护煤柱为研究对象。煤层埋深980 m，煤的容重为24 kN/m³，故上覆岩层的自重应力为24 MPa。

数值计算分为2 步。①在模型顶板施加恒定载荷，开展静载计算；②待静载计算平衡后，在模型顶板施加速度边界条件以模拟动载。

静载计算时，模型顶部设置为应力边界，上覆岩层自重采用均布荷载代替，底部为固定边界，顶部垂直应力取上覆岩层的自重应力(24 MPa)。在静载条件下，模型达到平衡状态，再在动载条件下开展计算。动力学计算时，数值模型顶部边界自由，底部设置为固定边界。

为了模拟动载扰动，假设上覆顶板岩层断裂产生含有多种子波的压缩波，而任何复杂的波形都可以由正弦波、余弦波或三角波合成。为简化动载扰动，选择正弦应力波作为输入波形。DUBINSKI[16]等现场监测发现应力波PPV 为0.03 m/s。而HADJIGEORGIOU 等[17]通过使用炸药诱发岩体振动所记录的应力波 PPV 不能大于3.3 m/s，因此，为了模拟不同震动强度的应力波，振幅设定为0.6、0.8、1.0、1.2、1.4 m/s。

2 动载扰动下煤柱失稳破坏变化规律

选取宽高比分别为4、6、8、10 四种模型，静载应力下平衡后，在煤柱模型的顶部分别施加振幅为0.6、0.8、1.0、1.2、1.4 m/s 的应力波（以应力波振幅表示动载强度），当模型再次达到平衡后，研究应力波输入保护煤柱破坏失稳规律。

2.1 煤柱尺寸对煤柱动力失稳的影响

由图3 可知，动载冲击后煤柱塑性区由表及里向深部扩展最终导致煤柱整体失稳。相比而言，随着煤柱尺寸增加，煤柱发生整体失稳所需的动载冲击强度具有一定的差异。以应力波振幅为0.6 m/s的模拟结果为例，以阐述此种差异。

图3 不同应力波振幅及不同尺寸保护煤柱塑性区分布规律Fig.3 Distribution law of plastic zone of protective coal pillars with different sizes and different stress wave amplitudes

煤柱宽高比为6 时，静载作用后，塑性区宽度为7～8 m，内部存在弹性核心。受到振幅为0.6 m/s 应力波冲击后，煤柱弹性核心完全破坏，塑性区贯穿整个煤柱，意味着煤柱整体失稳。相比而言，煤柱宽高比增加至8 与10 时，煤柱塑性区仅向深度扩展，内部依旧存在稳定的弹性核心，且弹性核心所占整个煤柱的比例由21%增加至47%。此外，对比发现，在相同静载作用下，不同尺寸的煤柱失稳启动所需的动载强度亦不同。模拟发现，随着煤柱宽高比由6增加至10 时，煤柱破坏失稳所需的动载冲击强度（以应力波振幅表征）由0.6 m/s 增加至1.2 m/s，提高了1 倍。由此说明，增加煤柱尺寸有助于提高煤柱的稳定性，减少动力失稳的几率。

2.2 动载冲击强度对煤柱动力失稳的影响

以煤柱宽高比为10 为例，分析动载冲击强度为0.6、0.8、1.0、1.2、1.4 m/s 时，煤柱失稳破坏特征。如图3c 所示，当应力波振幅从0.6 m/s 增加至1.0 m/s，煤柱塑性区由浅部逐步向着中部扩展，范围由8 m增加至13.5 m，煤柱内部弹性核区逐步减少。当动载冲击强度增加至1.2 m/s 时，煤柱整体产生塑性破坏，意味着煤柱整体失稳。由此说明，动载冲击强度越大，煤柱发生失稳破坏的概率也相对提高。

通过对比不同振幅应力波0.6、0.8、1.0、1.2、1.4 m/s 作用下，煤柱的失稳破坏特征，发现煤柱的稳定性与煤柱尺寸和动载冲击强度有关。相同条件下，尺寸越大，煤柱发生失稳破坏的几率越小；动载强度越大，煤柱破坏范围及程度越大，即保护煤柱失稳概率越大。

3 动载作用下煤柱突变失稳机制

3.1 煤柱承载力学模型

一般而言，煤柱载荷主要来源于上覆岩层重量和动载扰动。上覆岩层重量作为静载，使煤柱产生塑性区和弹性核区[18]，而后续动载冲击使塑性区向煤柱深处扩展，导致煤柱发生整体性破坏，内部弹性区减小而发生失稳，煤柱承载力学模型如图4 所示。

图4 煤柱承载力学模型Fig.4 Coal pillar bearing mechanical model

煤柱埋深为H，上覆岩层容重为γ，煤柱高度为Hm，煤柱宽度为B。假设外部动载q0均匀施加至煤柱，动荷载在某一时刻的运动状态可简化为静荷载[19]，动载简化系数可以表示为kd，则煤柱承受的载荷为：

煤柱塑性区与弹性区本构关系有所不同，弹性区煤体呈线弹性关系，塑性区内呈非线性关系，具有应变软化特征，煤柱塑性区内本构关系可以表示为

式中：E为煤柱的初始弹性模量；ε为应变；ε0为极限应变。

假设煤柱塑性区宽度为B0，则煤柱塑性区内载荷ps与煤柱压缩量u的关系可以表示为

式中：u0为煤柱最大压缩量，m。

弹性区内，煤体载荷pt与煤柱压缩量u的关系为

3.2 煤柱失稳的尖点突变模型

整个煤柱系统的总势能由煤柱自身的应变能和外部载荷所做的功组成。由式（4）可得煤柱塑性区应变能V1(u)为

由式（5）可得煤柱弹性区应变能V2(u)为

上覆岩层静载重力势能和外部动载荷做功的势能函数V3(u)为

根据上述力学模型的分析，动载作用下煤柱系统的总势能函数V(u)可以表示为

以u为状态变量进行突变理论分析，对V(u)求导得：

进一步简化，令：

将式（12）代入式（11）则有：

以x为状态参量，m、n为控制参量。由尖点突变理论可知，当m、n满足分叉集方程时，顶板-煤柱系统处于临界平衡状态。联立上式与其拐点方程3x2+m=0，消去状态参量x，得到分叉集方程：

由突变理论可知，当Δ＞0 时，系统处于稳定状态，煤柱不会失稳；当Δ=0 时，系统处于临界状态，受到微小扰动会发生失稳；当Δ＜0 时，是系统发生失稳的充要条件。将式（12）代入式（14）得：

煤柱失稳则方程（13）有解，必须满足Δ≤0，由此得到煤柱失稳的必要条件为：

3.3 煤柱最大压缩量计算

煤柱失稳判别式Δ中还含有一个与煤柱相关的未知量u0（煤柱最大压缩量），下面对煤柱最大压缩量u0进行求解。

煤柱的最大压缩量是由顶板的最大下沉量决定，如图5a 和5b 所示。考虑到煤柱的尺寸较大，拟采用弹性薄板理论研究顶板的最大下沉量[20-21]。图5c所示为直接顶-煤柱弹性薄板模型，其中，2a为煤柱长度，2b为煤柱宽度，q为直接顶上覆动静载等效均布载荷。直接顶可以视为四边固支的矩形板，其力学模型如图5d 所示。

图5 顶板-煤柱协同变形Fig.5 Theoretical process evolution diagram of elastic thin plate

该模型满足以下边界条件：

根据弹性力学中的伽辽金解法，将板的挠度展开为以下表达式：

式中：W为直接顶挠度；Ai为挠度系数；fi(x,y)为挠度函数，要求能满足薄板挠曲面的形状和几何边界条件，而且是连续的。

为便于计算，取一级近似函数：

显然该式能满足薄板挠曲面的形状和所有几何边界条件。根据伽辽金法对上式进行求解，得到挠度的表达式：

对挠度方程求极值，在板的中心处取到挠度的最大值：

煤柱最大压缩量等于直接顶最大挠度，如图5e 所示，即u0=Wmax，至此得到了煤柱失稳的最终判别式：

4 深井煤柱稳定性判别

4.1 深井煤柱稳定性判别方法

通过现场调研获得煤柱的物理力学参数，计算煤柱塑性区范围。

式中：k为应力集中系数；c0为黏聚力；φ0为内摩擦角；A为侧压系数。

计算煤柱最大压缩量：

最后，利用式（21）计算Δ值，进行稳定性判别。

4.2 工程背景

新河煤矿位于山东省济宁市唐口勘探区西南部边缘，目前开采水平为-980 m 延伸水平，延伸水平位于矿井的东南部，东西宽2 km，南北长2.4 km。新河煤矿正在开采6302 工作面，集中大巷分布在6302工作面以南，大巷之间的煤柱对大巷的稳定起到支撑保护作用。利用前文建立的煤柱稳定性判别方法，对新河煤矿6302 工作面南部轨道集中巷，6301 运输巷和6302 轨道巷之间的煤柱进行稳定性进行判别，如图6 所示。

图6 煤柱位置示意Fig.6 Position diagram of coal pillar

1) 计算煤柱塑性区范围：煤柱高2.9 m，黏聚力1.7 MPa，内摩擦角29°，埋深980 m，容重24 kN/m³，侧压系数1.12，应力集中系数1.97。将煤柱参数带入煤柱塑性区范围计算公式：

2) 计算煤柱最大压缩量：大巷保护煤柱长度46 m，宽度32 m，等效荷载3.75×107N，泥岩弹性模量6.5 GPa，泊松比0.4，顶板厚度6.63 m。

3) 利用尖点突变模型计算Δ值，并进行煤柱稳定性判别：

根据煤柱稳定性判别公式可知，Δ＝1.58＞0，该煤柱能保持稳定。

5 结 论

1) 数值模拟研究表明，煤柱受到外部动载扰动后，塑性区向中部扩展，弹性核心占比减小，导致煤柱破坏失稳。煤柱动力失稳与煤柱尺寸和外界动载冲击强度有关。增大煤柱的宽度能够提高煤柱的稳定性，减小煤柱破坏失稳的概率。相同条件下，煤柱破坏失稳的几率随着动载强度的增大而增大。

2）基于弹性薄板理论，建立了直接顶-煤柱弹性薄板模型，推导了直接顶的挠度方程，获得了煤柱的最大压缩量表达式。基于尖点突变理论，建立了煤柱尖点突变模型，推导了主要受塑性区宽度和煤柱压缩量影响的煤柱突变失稳判别式，为煤柱突变失稳提供了力学判据。

3) 基于煤柱突变失稳判别式，对新河煤矿6302工作面南部轨道集中巷，6301 运输巷和6302 轨道巷之间的保护煤柱进行了验证。结果表明，新河煤矿6302 工作面南部集中大巷保护煤柱能够保持稳定，与工程现场结果一致。