基于3DMine-Rhino-FLAC3D的多采空区对竖井的稳定性影响分析

郝勇浙，张 飞，王 昊

(内蒙古科技大学矿业与煤炭学院，内蒙古 包头 014010)

对矿产资源进行开采后，将在地下形成大量连续或非连续的采空区，采空区的存在破坏了岩体的原始应力平衡，使围岩应力重新分布。如果没有得到有效治理[1]，很可能造成采空区失稳，甚至引发岩体大面积的移动和冒落，采空区的失稳势必会对处在移动带内的竖井造成影响，对矿山的安全生产造成严重的威胁。伴随着数值模拟技术的不断进步，许多学者采用数值模拟技术对采空区及竖井的稳定性进行分析，取得了丰硕的成果。罗周全等[2]利用CMS探测及Midas-FLAC3D数据耦合技术极大提高了数值模拟的可靠性，为地下复杂采空区稳定性数值模拟分析开辟了新的途径；黄敏等[3]利用ANSYS和FLAC3D建立三维模型，对采空区的稳定性进行了分析，揭示了采空区周围的力学效应；王初步等[4]利用3DMine-FLAC3D耦合技术对采空区群进行研究，发现采空区应力叠加效应会使采空区发生破坏；赵永等[5]利用Mathews稳定图法计算出采空区的稳定数和等概率图，评价采空区的稳定性以及稳定概率。

本文以内蒙古某金矿为例，该矿开采深度接近1 000 m，形成了大量的采空区，竖井在移动带内且布置在矿体上盘，安全风险高，且对保安矿柱进行不规范回采，一旦采空区发生失稳，势必会对竖井造成威胁，形成重大安全隐患。 本文利用3DMine-Rhino-FLAC3D耦合建模，针对竖井的稳定状态进行分析，提出相应的治理措施并加以验证，由于盲竖井位于矿体下盘，相对较为稳定，此次不做分析。

1 工程概况

内蒙古某金矿13号矿体分布在矿区的东南部，是矿区的主要矿体之一。控制矿体长2 500 m，延深1 030 m，开采深度930 m(188标高)，平均厚度2.48 m。矿体呈大脉状-似层状产出，近东西向分布，矿体平均倾向183°，平均倾角55°，深部有变缓趋势。采用竖井加盲竖井联合开拓方式，明确竖井布置在矿体上盘且在移动带内，采用浅孔留矿法和全面采矿法进行开采。

截至目前，818 m中段以上已全部完成开采，818 m中段～188 m中段仍有矿房进行开采，由于大多数矿房在开采后没有及时得到处理，一些连续的采空区遗留在矿体的上部。采空区内留有顶底柱及间柱，现为矿区内最大规模的采空区域，拥有300多个小型采空区，受到工程扰动，采空区随时可能会塌落，对竖井造成严重威胁。竖井发生变形破坏会对运输、通风及安全生产造成重大影响，因此，对竖井周围的围岩进行研究分析是非常必要的，通过一定的治理措施改善围岩应力分布状态，降低应力集中，防止岩体过度变形，控制岩体运动幅度，可以消除安全隐患，避免灾害，这对矿山安全生产具有重要意义。

2 数值模拟

2.1 数值建模

矿山三维数字模型是三维数值模拟必要的数据基础，可以实现快速三维CAD向三维CAE的数据流通[6]。FLAC3D被广泛应用于岩土工程中，它能够进行岩石、土质和其他材料在达到屈服极限后经历塑性变形的三维空间行为分析，为求解三维问题提供了一种理想的分析工具[7]，但其建模能力较弱，尤其在建立复杂模型时工作量巨大。为了更加真实有效地构建矿体模型，基于三维网格堆砌法[8]并利用各软件的建模优势结合3DMine-Rhino-FLAC3D耦合建模，可解决复杂矿体、采空区的模型建立和网格划分等问题，使三维模拟与实际情况相对应，使模拟结果更加真实可靠。①使用3DMine软件和CAD软件的信息交互功能，将回采边界、矿岩分界面等信息导入3DMine中[9]，并矢量化线条，生成矿体，再利用生成DTM功能[10]、地质地形图的等高线，形成较为完整的三维地表模型。②利用Rhino及Kubrix模块[11]对建立好的矿体和地表进行x方向、y方向、z方向的边界圈定，依据主竖井及盲竖井的坐标及各中段间的位置关系，绘制5 m直径的井筒及半圆拱形巷道，并利用Griddle模块输出FLAC3D网格。③利用FLAC3D6.0的Model模块对围岩、矿体、井筒及采场进行分组，得到数值模拟的三维模型。

该金矿13号矿脉三维数值计算模型如图1和图2所示，该模型沿x方向长度为800 m，沿y方向长度为800 m，沿z方向的高度为1 200 m左右，由384 816个节点、676 921个单元组成，且沿z方向分为1 018 m中段～188 m中段，共20个中段。沿y方向上将各中段矿体划分为一系列的采场。分别对x方向(左右边界)、y方向(前后边界)及z方向的下部设置位移约束，z方向上部设为自由边界。

图2 矿体计算模型Fig.2 Ore body calculation model

2.2 本构关系

采空区和井筒的稳定性分析涉及围岩和周围矿体。 这些介质是弹塑性材料，适用于Mohr-Coulomb破坏准则的分析和研究。 因此，使用Mohr-Coulomb准则确定开挖后的矿体、采空区及竖井周围的应力，得到位移和塑性区的分布规律见式(1)和式(2)。

(1)

(2)

式中：σ1和σ3分别为最大主应力和最小主应力；C和φ分别为材料黏附力和内摩擦角；fs为破坏判断系数。当fs≥0时，材料处于塑性流动状态；当fs≤0时，材料处于弹性变形阶段[12]。在拉伸应力状态下，如果拉伸应力超过材料的拉伸强度，则材料将遭受拉伸破坏[13]。

2.3 岩石力学参数

根据该金矿所提供的岩石力学参数见表1。

表1 岩石力学参数Table 1 Rock mechanical parameters

2.4 计算方案

为准确模拟该金矿开采过程中竖井及其围岩的变化规律，根据开采深度和时间顺序，将开挖划分为22个步骤。 开挖顺序分别为主竖井及巷道、1 018 m中段、968 m中段、918 m中段、868 m中段、818 m中段、778 m中段、738 m中段、698 m中段、658 m中段、盲竖井及巷道、618 m中段、578 m中段、538 m中段、498 m中段、458 m中段、418 m中段、368 m中段、318 m中段、268 m中段、218 m中段、188 m中段。为简化计算模型，此次模拟不考虑竖井及巷道支护，研究不同开采深度下井筒和围岩内部力学状态的变化，包括岩体位移、塑性区和应力分布等。

3 采空区对竖井稳定性的影响分析

3.1 应力分析

图3为主应力变化曲线；图4和图5分别为垂直矿体走向上的围岩最大主应力和最小主应力的云图，随着矿体的逐渐开采，采场周围的原岩应力状态被打破，应力释放并进行重新分布。空区在上下盘岩层、顶底板形成次生应力场。由图3可知，随着矿体逐渐向下开采，各顶板的最大主应力和最小主应力均出现增大的趋势，开挖到188 m中段，最小压应力值为12.1 MPa，而最大压应力值为42.2 MPa。在顶板拐角处存在拉应力，远离顶板中央拉应力会逐渐减小。

图3 主应力变化曲线图Fig.3 Variation curve of principal stress

从图4和图5可以看出，随着向188 m中段矿体的逐渐开采，虽然采空区周围的围岩最大主应力和最小主应力均呈现增大的趋势，应力集中程度更大，但主竖井附近的应力状态仍然呈现出原岩应力状态，表明188 m中段～1 018 m中段之间矿体的开采对主竖井周围岩体的应力状态影响较小。

图4 采场及竖井最大主应力云图Fig.4 Maximum principal stress nephogram of stope and shaft

图5 采场及竖井最小主应力云图Fig.5 Minimum principal stress nephogram of stope and shaft

3.2 位移分析

图6为188 m中段～1 018 m中段之间的矿体开采完成后井筒附近的位移云图，随着矿体的逐渐向下中段进行开采，对井筒的影响范围是逐渐增大的。在采空区上方，围岩的变形呈现出拱形分布，距离采空区的距离越远，围岩的变形越小。

图6 采场及竖井位移云图Fig.6 Displacement nephogram of stope and shaft

图7为第1步到第20步开采过程中明竖井的9个监测点位移随着开挖步骤的变化规律，在188 m中段～1 018 m中段开采结束后，明竖井的最大变形量约为3.5 cm，最大变形位于658 m中段。目前，变形量在可控范围内，说明井筒处于安全状态，但主竖井的变形速率逐渐增大，由于岩石蠕变的作用以及爆破震动、地下水和节理裂隙等影响，若不及时对采空区进行治理，一旦采空区群发生失稳，可能会导致主竖井失稳，将会对矿山造成严重损失。

图7 采场开采不同开挖步骤下井筒变形曲线Fig.7 Shaft deformation curve under differentexcavation steps in stope mining

3.3 塑性区分析

图8为188 m中段～1 018 m中段之间的采场开采完成后竖井附近围岩的塑性区分布图，开挖的采场上盘围岩局部范围出现了剪切状态或拉伸状态，而主竖井附近的围岩仍然处于原岩应力状态，未出现剪切或拉伸状态，表明主竖井目前处于相对安全的状态。

图8 采场及竖井塑性区分布图Fig.8 Distribution of stope and shaft plastic zone

综上所述，在采场的每个中段进行开采后，一定范围内的围岩将发生拉伸或剪切破坏。其中，采掘采场中部浅段后，拉伸破坏主要发生在采空区周围；深部中段开采后，围岩破坏主要为剪切破坏；随着深部中段采场的开采，上部中段已经形成的塑性破坏区存在扩张现象，出现了一定的剪切破坏区，但尚未影响主竖井周围的应力状态，所以，主竖井处境较为安全。

4 治理措施

内蒙古某金矿818 m中段以上采场已经全部开采完毕，地表允许一定程度的冒落，这部分区域的采空区已经采用隔离封堵与自然冒落联合处理的方案进行治理。由于此次分析的采空区大部分位于保安矿柱之内，所以采用胶结充填对采空区进行处理，以保证竖井的稳定。由于采空区数量过于庞大，全充填的成本过高，考虑使用隔层充填的治理方案来节省成本。2种充填方案均采用向上充填的方式，全充填从188 m中段分16步充填至818 m中段，隔层充填则是从218 m中段分8步充填至818 m中段。

图9为采空区在两种不同治理方案处理后的井筒附近围岩位移分布云图，井筒附近的围岩变形情况基本一致，靠近采空区附近的围岩变形较大，越远离采空区，其变形越小；靠近采空区附近的围岩，其位移分布基本呈现拱形分布规律，采空区治理后，其分布规律基本与未处理前保持一致，采空区的顶板最大变形值约为7 cm，表明治理后采空区顶板的变形速率较小，充填体能有效支撑上盘围岩，保证采空区的稳定性，进而保证主竖井附近围岩的稳定性。

图9 充填后采场及竖井位移云图Fig.9 Displacement nephogram of stope and shaft after filling

图10为2种治理方案在采空区处理前和处理后主竖井各个监测点的变形情况，其中，第1步～第20步为采空区未处理情况下的各监测点变形情况，20步以后为采空区在不同处理阶段各监测点的变形情况。采空区处理前，主竖井各监测点的变形值较大，变形速率也较大，最大变形值约为3.5 cm，2种方案处理后的主竖井的变形速率均呈现减小的趋势，表明采空区处理后井筒的变形均得到一定的控制，说明2种充填方案均能对上盘围岩起到支撑作用，其稳定性可得到一定的改善。但隔层充填方案的变形速率比全充填方案的变形速率要高，变形量也较高，最大变形值为3.74 cm，比全充填高0.04 cm。所以在条件允许的情况下，建议使用全充填的治理方案，与监控系统配合，实现采空区管理过程中地压危险区和危险程度的实时监测、预警和预报，可为矿井下危险源区域和危险程度分析提供可靠数据支持，保障井下作业人员和设备的安全。

图10 不同开挖与充填步骤下主竖井变形曲线Fig.10 Deformation curve of main shaft under different excavation and filling steps

5 结 论

1) 通过对该金矿的采空区及竖井进行调查，以三维网格堆砌法为基础，提出了3DMine-Rhino-FLAC3D耦合建模的建模思路，能够更加快速、准确地构建采空区及竖井的三维模型，并就采空区对竖井稳定性的影响进行模拟分析，发现在1 018 m中段向188 m中段开采过程中虽然在采空区周围的围岩应力有增大的趋势，但在开采过程中对竖井周围应力影响较小。

2) 通过对位移和塑性区的计算结果可知，虽然空区周边发生拉伸破坏和剪切破坏，且采空区周边塑性破坏区存在扩张的现象，但主竖井附近的围岩仍处于原岩应力状态，并未出现塑性区。在开挖过程中，主竖井各中段发生了不同大小的变形，最大变形量约为3.5 cm，位于658 m中段，虽然变形在可控范围内，但变形速率不断增大，竖井仍有可能发生开裂破坏。

3) 通过对隔层充填和全部充填2种治理方案进行模拟计算，发现2种治理方案在充填后的竖井变形速率均明显减小，趋于平缓，说明充填体能够有效支撑上盘围岩，使井筒变形的到控制，虽然全进行胶结充填成本过高，但隔层充填的变形速率与变形量大于全充填的治理方案，所以推荐使用全部充填的治理方案。