深部崩落法开采地表沉降特征数值模拟研究

王洛锋

(洛阳栾川钼业集团股份有限公司，河南 洛阳 471500)

0 引言

崩落法开采引起的地表沉降是一个重大的矿山安全和环境问题，研究矿体开采后岩层的移动规律，对矿山安全高效开采及社会经济有着重大意义。对于采用崩落法开采的矿山，矿体开采后，原岩的应力平衡会被打破，将会引起围岩发生位移、开裂和冒落等。由此产生的地表错动、塌陷将会影响到工业场地等地面设施的布置，同时也会引起生态环境的损坏以及征地等问题。因此，能够有效的预测崩落法开采引起的地表塌陷范围，对矿山规划、风险管理与评估非常重要。许多学者针对地表沉降方面已开展了大量研究。如王永嘉应用离散元法和边界元法对开采沉陷中的应力、变形进行了研究[1]；唐春安给出了线弹性有限元法[2]；张东明、尹光志等人探讨了煤矿开采中岩层移动的分形特征[3]；郭增长建立了极不充分开采条件下地表移动的预计方法[4]；黄平路、陈从新采用多种数值模拟手段综合分析方法，研究了露天和地下同时开采引起的岩体移动规律[5]，均取得了一定的成果。

前人的研究内容多集中于埋深比较浅的情况，随着开采深度的增加，深部岩体变形特征将发生转变[6]，深部开采引起的岩层移动规律与浅层开采不同。因此，研究深部矿体开采时地表沉降规律具有重要的现实意义。基于岩土力学理论和有限差分软件FLAC3D，对深部矿山采用崩落法开采引起的地表沉降特征进行研究。

1 工程背景

某铁矿地表标高为0～30 m，矿体为急倾斜极厚矿体与厚矿体，主矿体走向270°～310°，倾角70°～80°，倾向西南，矿体厚度80～270 m，平均厚度为160 m，矿体延伸较大，深处未见尖灭。采矿方法为无底柱分段崩落法，最低开采标高为-950 m。矿体回采顺序在矿体走向上由厚大部位两端推进。根据上述工程背景及矿山产量要求，设计开采方案概念图为图1，退采的过程中保持4～5个水平同时开采，开采分段保持相对平缓的“V”形向下逐分层推进，共分17步采完。

图1 开采过程概念图

矿床基岩属于坚硬的、次硬的层状、块状工程地质岩组，但岩组结构较复杂，各类结构面较发育。浪子山组底砾岩为其直接顶板，岩石较完整，但其厚度不均，一般为2.72～50.42 m。硅化大理岩、绿泥绢云石英片为其间接顶板，岩石硬脆，节理裂隙较发育，钻孔揭露岩心较破碎，岩体稳固性较差。铁矿层下盘围岩绿泥片岩，岩石较软、破碎，节理裂隙较为发育，多见破碎带。矿区工程地质条件属中等类型。

2 计算模型与计算过程

2.1 数值模型的建立

该矿山矿岩种类较多、性质复杂，绿泥石英岩与绿泥石英片岩构成矿体的直接顶、底板，根据地质剖面图与采矿方法参数，建立的计算模型为图2～3。模型尺寸为2.8 km(x)×3.8 km(y)×2 km(z)，在开采范围区域及其附近采用较为精细的网格，边界区域采用相对稀疏的网格，模型共划分为1 742 626个单元，292 222个节点。模型的底面采用全约束，四个侧面采用法向约束。断层采用接触面(INTERFACE)模拟。

图2 计算模型

图3 岩层分布图

2.2 本构模型的选择及力学参数的确定

在此次模拟计算中，远离扰动区的位置采用摩尔-库伦本构模型，靠近扰动区的位置采用岩体劣化本构模型(RDM)[7]，该模型在摩尔-库伦的基础上优化提出，该本构模型是基于破坏范围内的强度参数随着岩体破坏程度的不同而发生相应劣化的基本特征，认为岩体破坏范围内的岩体力学参数是等效塑性应变的函数，如式(1)所示：

(1)

强度参数的确定采用了基于Hoek-Brown强度参数的换算方法，基于岩体质量分级结果RMR，根据经验公式换算成GSI，采用Hoek-Brown强度准则计算获得岩体力学参数，见表1。

表1 矿岩力学参数值汇总表

2.3 计算方法及步骤

此次数值模拟采用国际通用的大型三维岩土计算软件FLAC3D，该方法是基于Cundall P.A.博士提出的一种显式有限差分计算方法。经多年发展，该软件已经成为岩土工程领域的通用软件。该软件计算过程中不需利用迭代来满足相应的本构关系，而是根据应变与应力的函数关系，使应力随着应变的变化而变化，因此在处理复杂岩体条件下的岩石工程开挖计算具有独特的优势。

计算初始应力场按构造应力场考虑，根据地应力测量结果，垂直方向地应力与自重应力γH基本相同，水平构造应力σx=λ1σz，σy=λ2σz，λ1、λ2为侧压力系数，根据地应力测试结果，分别取值为1.5、1.7。按照图1的计算方案进行计算，计算矿体崩落后后围岩应力场变化及相应的围岩变形，直至完成17个开采步，即从-740 m到-950 m水平的开采。计算步骤如下：①赋于模型相应的岩石力学参数及边界条件，模型的侧面施加法向约束边界条件，模型的底部设置x、y、z三向固定边界条件。②在重力场作用下计算初始地应力场，计算平衡后对计算得到的位移、塑性区进行清零，模拟开采前的岩体力学状态。③基于选取的本构模型，设置与崩落模拟相关的计算参数。④根据开采过程，以每个开采分段作为一个计算步进行计算，获得每个开采分段的应力场、位移场及塑性区，直道计算趋于平衡。⑤对上述开采步骤进行重复计算，待整个矿体开采结束后停止计算。

3 数值模拟结果分析

从变形特征来看(见图4～6)，崩落区形态的基本特征是崩落边界呈陡立状，意味着地表沉陷区范围相对较小，下盘断层的存在对变形起到一定的“隔断”作用。在第3步(-776 m)及其之前的开采过程中，崩落区形态总体呈下宽上窄的“烟囱”形。随着开采的继续，如第4～10步开采过程中，崩落区轮廓总体保持陡立，从第12步(-884 m)开始，开采区上下盘一带开始出现一定程度的横向变形或崩落，并在后期阶段对地表变形和沉陷范围形成一定的影响。随着开采范围的增大，崩落区形态会发生一定变化，变形逐渐增大并影响到地表的沉陷和变形分布范围。图7为开采过程中地表沉陷曲线，可以看出，下盘岩体位移发展到一定程度后，位移变化速度将变缓慢，上盘岩体随着矿体的继续开采，岩体移动范围有变大的趋势，总体而言，受矿体赋存产状等因素影响，矿体上盘对地表的影响范围要大于矿体下盘对地表的影响，开采结束后地表最大位移约为14 cm。

图4 第3步开采完成后围岩变形

图5 第12步开采完成后围岩变形

图6 第17步开采完成后围岩变形

图7 地表沉降曲线

从应力特征看(见图8～9)，矿体开挖后，采场上部在开采前期会出现一定范围的“应力拱”，随着矿体开采的进行，应力承压拱的范围逐渐扩大，承压拱可以对上覆岩体起到一定的支撑作用。但由于应力分布会随着开采的进行而不断的发生变化，最后将会引起应力承压拱的力学条件发生破坏，引起上覆围岩岩破坏，地表发生塌陷。

图8 第6步开采完成后最大主应力云图

4 结语

根据崩落法开采过程，基于FLCA3D数值模拟软件，获取了深部矿体崩落法开采位移及变形规律，研究成果具备一定的工程参考价值。主要结论如下：

1)深部矿体采用崩落法开采后，对于急倾斜厚大矿体而言，沉陷类型基本表现为筒状陷落，断层的存在会对位移的传递起到一定“隔断”作用，移动角要比浅部传统方法确定的大。

2)矿体开采后，采场上方存在由高应力组成的“复合应力承压拱”，其位置随崩落的进行而不断演化，承压拱的存在会对上覆岩层起到一定支持作用，最终成拱机制破坏，导致地表的沉陷。