考虑围岩应变软化及采空区接触的深部煤层底板破坏分析

尹尚先，孟浩鹏，钱双彬

(1.华北科技学院 安全工程学院，北京 101601；2.华北科技学院 理学院，北京 101601)

随着时代发展，数值方法在岩土工程领域的应用逐渐普及，其中FLAC3D有限差分数值方法优点众多，在采矿领域得到广泛认可与应用，众多学者借助FLAC3D围绕煤炭安全开采进行了大量研究。董书宁等[1]在改造奥灰顶部岩层段的判别准则研究中，利用数值计算分析了采深、采高和采宽等尺寸效应对底板破坏带的影响；刘伟韬等[2]通过正交试验设计，进行了7 个主控因素对底板破坏深度的影响研究，并对主控因素进行排序；刘新民[3]采用现场和模拟结合的方法进行沿空留巷对底板破坏深度的研究，认为无煤柱式的沿空留巷开采技术不会对底板破坏深度造成较大影响；朱斯陶[4]、朱广安[5]、田雨桐[6]等的研究表明，数值模拟是研究采动影响下断层活化规律的重要手段之一，可实现采动影响下对断层多方位的定量分析；朱庆伟[7]、甘智慧[8]等运用数值方法对采动影响下覆岩结构演化和地面沉降进行了研究。上述研究成果在数值计算中一般采用摩尔-库伦本构模型，该本构模型为理想弹塑性模型，对围岩塑性屈服后的状态无法准确描述，并且对拉格朗日法遵循连续介质假设而导致网格发生大变形但节点不接触的固有缺陷[9]未进行深入研究，这2 个问题可能拉大数值模拟同实际情况之间的差距。大变形条件下不接触的固有缺陷，导致无法模拟采空区顶板垮落后顶底板接触的应力传递现象，即“采空区不接触”现象。观察模拟工作面回采结果发现：顶底板接触后并不进行接触计算；顶底板应力基本处于泄压状态。真实情况是采空区在周期来压过程中，基本顶发生周期性垮落与底板接触，采空区的应力恢复随时间变化呈指数函数关系[10]。因此，利用FLAC3D研究煤层回采后底板破坏时，应对围岩应变软化和采空区接触进行考虑。

笔者将针对FLAC3D模拟工作面回采中本构模型的选择和采空区不接触的固有缺陷进行研究，以河北开平煤田林西矿2023 工作面底板导水裂隙带实测为工程背景，建立考虑应变软化和采空区接触的工作面回采模型，并结合力学分析对模拟结果进行解释，达到提高工作面回采数值计算准确性的目的。

1 工程背景及底板破坏深度实测

1.1 矿区概况

林西矿位于河北开平煤田东南翼，地质构造以褶皱为主，开平主向斜穿过井田东侧深部，西部有杜军庄背斜和黑鸭子向斜；井田内大型断层较少，小型断层较发育；井田内揭露地层由老到新为：奥陶系、石炭系、二叠系和第四系。石炭系、二叠系为含煤地层，煤系基底为奥陶系灰岩(简称奥灰)。林西矿目前主要生产水平为11 水平(−850 m)和12 水平(−1 000 m)。

12 煤位于二叠系下统赵各庄组，煤层底板至奥陶系自上而下为：石炭系上统开平组，上以K6 灰岩顶界面与二叠系下统赵各庄组分界，下以唐山组K3 灰岩顶界面与开平组分界；石炭系上统唐山组，上以K3 灰岩顶界面与开平组分界，下以G 层铝铁质泥岩的底界面与奥陶系石灰岩呈平行不整合接触。赵各庄组为主要含煤组，厚度33.55～61.20 m，平均48.35 m，含11、12 两层可采煤层；本组岩性顶部主要为黑色泥岩和灰色砂岩及褐灰色砂岩，其中砂岩向深部逐渐发展为黑色泥岩；中部及底部主要为灰色细−粗砂岩，浅部间有3～4 层砂砾岩，向深部砾岩直径逐渐变小，砂岩粒度也逐渐变细且大部分砂岩为泥岩所替代。开平组层厚55.09～91.04 m，平均76.81 m，主要岩性为黑色粉砂质泥岩，砂岩次之，其中砂岩比例由浅向深逐渐减少。唐山组层厚64.44～80.55 m，平均70.18 m；地层岩性除K1、K2、K3 灰岩和G 层铝铁质泥岩外，主要为黑色-深灰色泥岩和灰色砂岩。

根据林西矿深部ZK7 奥灰水位观测孔钻孔水位及12 煤层底板至奥灰顶界面间距计算得到：当12 煤层底板高程小于−936.3 m 时，工作面回采期间突水系数超过0.06 MPa/m。按照《煤矿防治水细则》，如果突水系数超限问题不解决，深部区域将无法进行安全带压回采，严重影响矿井采掘接替，威胁矿井生存。因此，有必要对12 煤层底板采动导水破坏带深度进行实测及底板破坏分析等工作。

1.2 工作面概况

林西矿深部2023 工作面开采12 煤层，位于林西井田杜军庄背斜构造块内，煤层走向变化较大(N11°EN36°E)，煤层厚度0.8～2.7 m，平均2.0 m；煤层倾角17°～21°，平均20°；走向长约661 m，倾斜长约93 m；地面高程30 m，开采高程−842.8～−884.0 m。直接顶为炭质泥岩，厚约3.13 m；基本顶为泥岩，厚约2.08 m；直接底为泥岩，厚约0.7 m，老底为粉砂岩，厚0.5～4.6 m。

1.3 底板破坏深度实测

为观测林西矿2023 工作面底板采动导水破坏带发育情况，在2023 工作面西侧2023-2 工作面回风巷设计D01-1、D01-2 钻孔，钻孔钻至2023 工作面法向向下30 m 左右位置，如图1 所示，随后进行分段压水试验，记录不同位置流量稳定后的漏失量，绘制钻孔的漏失量曲线(图2)。观察钻孔的漏失量曲线发现，在垂距24 m 以后，漏失量下降明显，注水压力在漏失量稳定后回弹明显且接近初始注水压力，说明该段裂隙不发育且贯通性差。在压水试验基础上，为直观了解底板破坏情况以及对压水试验结果进行验证，对D01-2钻孔进行钻孔成像，获取成像数据如图3 所示，图中距离均已换算为距煤层底板垂深。

图1 2023 工作面底板破坏深度实测钻孔布置Fig.1 Layout of boreholes for the failure depth in the floor of working face 2023

图2 2 个钻孔漏失量曲线Fig.2 Two boreholes leakage curves

图3 D01-2 钻孔不同深度成像Fig.3 D01-2 borehole imaging at different depths

根据D01-2 钻孔成像结果显示，图3a 中14.2 m处岩性为灰色细砂岩，处于底板破坏带边缘，该处岩层完整性较好，但存细微层状裂隙；图3b、图3c 中15.9～19.3 m 处，岩性为细砂岩，这2 处岩层及其之间的岩石破碎严重，纵横裂隙发育明显；图3d 中26.6 m处，岩性为细砂岩，该处岩层及以下岩层完整，未见明显裂隙。认为D01-2 钻孔漏失量曲线在30 m 处的回弹是由于压水试验中的封堵装置压力过大致使原生裂隙张开导致。

综合D01-1、D01-2 孔压水试验实测结果和钻孔成像数据，并参考邻近赵各庄矿近似埋深的1237、2137 工作面正常底板实测采动导水破坏带深度23、25 m[11]，最终确定林西矿2023 工作面正常底板采动导水破坏带深度为24 m。

2 底板破坏力学分析

运用朗肯土压力理论对煤层底板破坏进行定性分析[12-13]，图4 表示具有半无限平面的煤层开采走向剖面图，将底板塑性破坏划分为Ⅰ区(主动区)、Ⅱ区(过渡区)、Ⅲ区(被动区)3 个区域。图4 表示回采过程中煤层底板受力状态的应力圆与底板强度包线之间的关系。

图4 煤层底板破坏分区Fig.4 Coal seam floor failure zones

在原岩应力状态下，煤层下深度为z处单元体的应力为竖向应力σv(σz)、最大水平应力σH、最小水平应力 σh，取σH=σ1(最大主应力)、σv=σ3(最小主应力)的情况，用图5 中的应力圆①表示，此状态下应力圆距强度包络线较远，底板处于弹性平衡。随工作面的推进，底板地应力受工程扰动，在超前应力作用下 σv将逐渐增大，假设σH保持不变，此时应力圆的半径先减小后增大，如图5 中方向向右的箭头所示，若剪应力达到底板抗剪强度，应力圆与强度包络线相切，该处底板达到被动极限平衡状态，如图5 中应力圆②所示，此时σH=σ3、σv=σv1=σ1(σv1为应力圆②状态下最大主应力)，根据摩尔−库伦理论可知，当工作面前方底板达到或超过极限平衡状态并随应力的持续作用，Ⅰ区发生塑性变形，伴随体积膨胀以压力形式通过Ⅱ区向采空区方向(Ⅲ区)传递，导致底鼓同时形成连续的滑移面。当工作面推过之前超前应力作用的底板位置后，最大、最小主应力发生转变，由于工作面的推进，σv减小，假设σH保持不变，此时应力圆的半径先减小后增大，如图5 中方向向左的箭头所示，直至应力圆与强度包络线相切，该处底板达到主动极限平衡状态，如图5中应力圆③所示，此时σH=σ1、σv=σv2=σ3(σv2为应力圆③状态下最小主应力)，根据摩尔-库伦理论可知，当工作面前方底板达到或超过极限平衡状态并随应力的持续作用，Ⅲ区发生塑性变形，导致采空区底板隆起，同时形成连续的滑移面。

图5 煤层底板极限平衡状态Fig.5 Limit equilibrium state of coal seam floor

为比较底板的被动极限平衡状态和主动极限平衡状态，根据极限应力圆与强度包络线之间所得的关系式[14]：

式中：φ为岩石内摩擦角；c为岩石纯剪切强度(黏聚力)。

将应力圆②和③的表达式代入式(1)可得底板被动和主动极限平衡状态之间的关系：

需注意：①实际煤层开采工作面斜长有限，同时工作面斜长与底板采动破坏带深度密切相关，导致对底板破坏进行定量分析存在困难，但采用半无限平面对底板破坏进行定性力学分析是可行的；② 对煤层下深度为z处单元体最大水平应力 σH保持不变的假设同实际不符，Ⅰ区的 σH应当随顶板周期来压而发生周期性变化，Ⅰ区在竖向应力作用下，应力必然向四周传递，由于采空区这一临空面的存在，Ⅱ区、Ⅲ区的 σH随Ⅰ区的 σH的变化而变化；③从回采过程中底板受力状态的应力圆与底板强度包络线之间的关系可知，达到应力圆②的状态比应力圆③要困难，即Ⅰ区达到极限平衡状态较难。

3 数值模拟论证

3.1 本构模型的选择

底板破坏模拟一般采用摩尔-库伦本构模型，此模型为理想状态的弹塑性模型，不考虑黏聚力、内摩擦角、抗拉强度等材料参数随塑性变形的变化情况[15]，同实际地质材料受力特征不符。此模型在底板塑性屈服前能较好地反映底板变形，但在底板发生塑性屈服阶段开始后，同实际破坏有较大差异。底板塑性屈服后，在应力作用下呈应变软化行为，底板在地应力的作用下产生微裂纹及岩体的相对滑动，底板强度将不断降低并且越来越缺乏弹性，直到破坏以及剪切带的产生。利用FLAC3D建立应变软化摩尔-库伦地层模型，能够对围岩塑性破坏后的力学状态更准确表述。

实际岩石峰后的应变软化过程的弹塑性刚度矩阵为一个不定矩阵[16]，导致应变软化问题求解困难。为避免这一情况，可将岩石峰后应变软化过程简化为一系列的脆塑性过程[17-19]。FLAC3D内置的应变软化模型为基于经典弹塑性理论将实际的应变软化过程中黏聚力、内摩擦角、剪胀角与塑性剪切应变的函数近似为一组首尾相连的分段线性函数的模型。

基于岩石软化相关文献[20-21]，将其实验数据拟合为岩石峰后黏聚力、峰后内摩擦角随塑性剪切应变εs的指数函数：

为验证拟合函数的合理性，将函数嵌入应变软化本构模型，对参考文献中的泥岩进行单轴压缩数值模拟。模拟采用单轴压缩试验常用直径(D)∶高(H)为1∶2 的圆柱体进行试验，为较好地模拟真实试验和呈现应力-应变全过程曲线，对试件端部施加恒定速度代替试件受压情况。根据上述条件进行2 种本构模型的试验，得到σ-ε曲线(图6)。屈服前，2 种材料的σ-ε曲线一致且基本符合线弹性；屈服后，摩尔-库伦材料与应变软化材料的σ-ε曲线明显不同。应变软化材料的σ-ε曲线同真实试验有较高的吻合度，可较好地描述岩石屈服软化后的力学特征，说明根据式(3)、式(4)建立的应变软化模型是合理可行的。

图6 2 种本构模型应力-应变曲线Fig.6 Stress-strain curves of two constitutive models

3.2 模拟过程

3.2.1 初始模型建立

数值模拟的工况条件以林西矿12 煤2023 工作面为背景，由于12 煤层厚度0.41～8.48 m，煤层含夹矸0～1 层，夹矸厚度0.10～0.31 m，结构较简单，埋深可至1 000 m 以下，故模拟工作面采高4 m，倾向与走向长度100 m×800 m，埋深1 000 m。地层信息参考林西矿深部ZK7 奥灰水位观测孔钻孔信息，对地层倾角、岩石力学参数相近及薄岩层进行适当简化，最终模型尺寸长×宽×高为1 000 m×300 m×240 m，剖分网格数量804 000 个，节点数量833 748 个(图7)，数值模型的岩石力学参数见表1。

图7 煤层底板破坏数值模型Fig.7 Numerical model of coal seam floor failure

表1 数值模型岩石力学参数Table 1 Rock mechanical parameters of the numerical model

本构模型选用本文提出的应变软化模型，对模型施加10 m/s2竖直向下的重力加速度，在模型顶部施加22.618 MPa 应力代替未建模的上覆岩层，底部采用竖直位移约束，四周采用水平位移约束，并在四周施加随深度增加的侧向应力。侧向应力的大小参考沉积岩的应力分布规律[22]，根据模型煤层埋深1 000 m，得垂直应力(σv)∶最大水平主应力(σH)∶最小水平主应力(σh)为1∶1.133∶0.758，为研究原岩应力对结果的影响，在模拟中施加2 种相反地应力分布，设置初始地应力分布情况：σZ∶σX∶σY=1∶1.133∶0.758、σZ∶σX∶σY=1∶0.758∶1.133(X、Y、Z为模型坐标方位)。

3.2.2 考虑采空区接触的方法

实际煤层开采中，基本顶随周期来压垮落，垮落后的顶板与底板接触，发生顶底板之间的应力传递。但在FLAC3D模拟工作面回采中发现，顶底板在应力作用下接触后，并不会进行应力接触计算，在大变形模式下，甚至可清楚观察到顶底板发生交叉的现象，这与实际情况不符，导致模拟与实际产生巨大偏差，所以在模拟工作面回采中需要对采空区进行接触模拟。为体现真实回采中应力变化，开挖步距参考真实工作面周期来压步距，模拟过程中工作面以20 m 为步距循环开挖。煤层回采后，回采区域由“应变软化”模型转为“空”模型。将实际采空区顶板垮落后堆积的碎石假设为弹性整体，利用“弹性(各向同性)”模型替换“空”模型，达到模拟采空区顶板垮落后顶底板应力接触的目的。在弹性体参数的确定上，由于采空区的碎石是顶板垮落产生，将其视为裂隙发育的弹性整体，其弹性模量将大幅衰减，泊松比有所上升[23]。

数值计算时，工作面回采3 个循环后，模拟顶板垮落后的顶底板接触，即“空”模型以20 m 为步距循环向“弹性”模型转化，为避免弹性体对侧向产生应力传递，在切眼、侧帮及终采线附近不改变“空”模型。在工作面中心顶底板位置分别布置测点，记录回采全过程竖直方向应力值和竖直位移量(图8、图9)。根据测点记录发现，模拟采空区顶板垮落后的接触与对采空区不做处理的结果具有显著差异，考虑采空区接触的结果更贴近实际。考虑采空区接触的采空区顶底板应力得到一定程度的恢复，而不考虑接触的采空区顶底板应力基本处于泄压状态；考虑采空区接触的采空区顶底板竖直位移量明显小于不考虑接触的情况，并且考虑采空区接触后的底板位移量表现出小幅的回落。

图8 煤层顶底板竖直方向应力Fig.8 Stress in vertical direction of coal seam roof and floor

图9 煤层顶底板竖直位移Fig.9 Vertical displacement of coal seam roof and floor

采用“应变软化-空-弹性”模型转变的方法，达到模拟采空区顶板垮落后应力传递的效果，弥补了以往煤层开采模拟中采空区顶底板不接触的固有缺陷。

3.3 模拟结果及分析

运用自定义应变软化本构关系和考虑采空区接触的数值方法进行目标工作面的回采模拟，对地应力σZ∶σX∶σY=1∶0.758∶1.133 模拟过程中(回采40、100、800 m 平衡后)的塑性区(图10)和累计塑性剪切应变率大于0.01 的区域(图11)进行切片展示。

图10 不同回采距离塑性区分布Fig.10 Plastic zone distribution in different mining distance

图11 不同回采距离塑性剪切应变突出区Fig.11 Strain-shear-plastic outburst zone in different mining distance

1) 底板塑性区分析

模拟回采0～100 m 过程中，塑性区深度迅速增大，后随回采的进行，塑性区深度基本稳定在23 m 左右，同实际底板导水裂隙带深度的实测结果一致；底板塑性区上部的状态为“shear-p、tension-p”，为剪切屈服和张拉屈服共存状态，分布形态随回采呈周期性分布；底板塑性区下部的状态基本为“shear-p”，为剪切屈服状态，将上部“shear-p、tension-p”状态包围。通过结合该区域应力及位移分布，底板塑性区上部的最小主应力呈拉应力，下部呈压应力；在底板2 种塑性状态交界附近的位移量有一定突变。结果表明，底板破坏区域上部为剪切、拉张交互破坏，下部为剪切破坏。

2) 塑性剪切区分析

底板塑性区全区包含塑性剪切状态(图10)，选取累计塑性剪切应变率大于0.01 的区域进行显示(图11)，其分布形态为斜向采空区的半包围面状结构。

截取2 种极限地应力条件下的3 个开采循环距离进行底板破坏分析，对塑性区、塑性剪切应变率大于0.01 进行整合处理(图12)。发现原岩应力的改变对底板破坏规律几乎没有影响，塑性剪切应变突出区域将剪切、拉张交互破坏区(“shear-p、tension-p”)同剪切破坏区(“shear-p”)划分开，即剪切带内侧为剪切、拉伸交互破坏，外侧为剪切破坏。一般认为，岩石沿最大有效剪应力面形成剪切破裂面，而塑性剪切应变集中区域与最大有效剪应力集中区域基本一致，故认为滑移面即剪切破坏面是沿塑性剪切应变集中区域分布。通过利用优化后的数值方法得到的结果可知，煤层回采后底板破坏类型可分为剪切和拉张交互破坏、剪切破坏2 种类型，并可根据塑性剪切应变集中度，对底板滑移线进行三维可视化显示，底板滑移面呈斜向采空区的半包围面状阵列分布。

图12 2 种地应力下局部模拟结果Fig.12 Local simulation results under two ground stresses

将数值结果与本文提到的底板破坏力学分析结合，工作面回采前，底板岩体处于地应力平衡状态，受采动影响底板初始应力状态被打破，在采空区前方底板(Ⅰ区)产生超前支撑压力，由于采空区这一临空面的存在，Ⅰ区的大部分应力和位移通过Ⅱ区向Ⅲ区传递，致使采空区底板整体处于高围压、低轴压状态，底板破坏形式整体表现为塑性剪切破坏。对塑性剪切应变率较高区域进行显示，该区域将底板塑性破坏区分为上下部分。朗肯土压力理论中3 个破坏区是根据塑性剪切程度进行划分，认为模拟结果中底板塑性区可根据塑性剪切应变率的大小进行划分，划分结果可与朗肯土压力理论中3 个区进行对应。本文底板塑性区中塑性剪切应变率突出的区域(滑移面)为Ⅱ区和Ⅲ区之间的分界，Ⅰ区未显现，同前文“Ⅰ区达到极限平衡状态较难”对应。Ⅱ区的破坏形式为剪切破坏，Ⅲ区的破坏形式为拉张和剪切的交互破坏。

4 结 论

a.运用朗肯土压力理论并结合考虑围岩应变软化和采空区接触的数值方法研究了河北林西矿深部底板破坏特征，根据塑性剪切应变率的变化，对底板滑移面实现了三维显示，并将底板塑性区与朗肯土压力中的主动区、过渡区和被动区相对应，其中过渡区、被动区破坏形式分别为剪切破坏、拉张与剪切的交互破坏。

b.提出的考虑围岩应变软化和采空区接触的FLAC3D数值方法，对煤层开采数值模拟实现了优化。在本构模型的选择上，建立更贴合实际的应变软化本构关系；对以往模拟中采空区顶底板不接触的固有缺陷，采用“应变软化-空-弹性”模型转变的方法得到解决。该方法为煤层开采及需要考虑开挖后接触的大变形工程的数值计算提供一种更贴合实际的模拟思路。

c.应变软化本构模型能够让计算结果更贴合实际，但该模型需以大量岩石峰后黏聚力、内摩擦角衰减的数据为基础进行建立，因此，还需对不同岩石峰后的力学现象进行深入研究。