摘要：在对五道羊岔矿区岩体宏观结构调查及室内岩石力学试验的基础上，采用FLAC有限差分法对矿区开采诱发的地表沉陷变形进行了数值模拟研究，得出矿区5个采区在不同采矿技术参数下的围岩移动及地表变形的规律。基于数值模拟结果，确定了5个采区在开采完成后的地表沉陷边界，并估算了各采区的地表塌陷体积。
关键词：采矿工程,地表变形,FLAC数值模拟,沉陷边界,塌陷体积
在矿山生产中，人们往往希望最大限度地对矿产资源进行开采;另一方面又要满足矿山在开采期间中的稳定性。然而在工程实践中，矿产资源的大规模、高强度开采，已对生态环境造成了严重的负面影响，导致地表下沉、塌陷等，严重影响矿山及地表构筑物的稳定性。因此，对不同采矿方法、采矿工艺及技术参数下矿区地表沉陷、变形的预测具有重大的研究意义和实践指导意义[1]。
地下开采诱发的岩层移动及地表沉陷是一个复杂的物理、力学变化过程[2]。国内外学者提出了许多关于开采沉陷的预测方法和沉降模型，如典型曲线法、剖面函数法、影响函数法、下沉网格法等[3]。20世纪80年代，伴随着随机介质理论的引入，我国学者在此基础上建立了地表移动预计的经典概率积分法，是我国较为成熟和应用最为广泛的预测方法之一[4-5]。然而，概率积分法在对岩层移动力学机理上的解释存在不足，且关键参数须结合实地观测数据确定，这给开采沉陷的预计工作造成了一定的难度[6]。
随着离散元、有限元、有限差分及边界元等各种算法和程序的问世，应用数值分析方法研究开挖引起的地表扰动变形得到了长足的发展。在各类数值计算方法中，有限差分法应用最为广泛。它可模拟不同采矿工艺[7-8]和采矿方法[9-10]条件下地表的沉陷变形，从而优化观测站布置[11]，圈定地表移动带范围[12]，进而为采动灾害防治及地表建(构)筑物的保护措施提出合理建议。
以五道羊岔矿为例，采用数值计算方法研究了地下采动作用下矿区地表变形规律，圈定了地表沉陷带的范围，并估算开采诱发的塌陷体积，分析评价其对地表环境的损害，为该方法在地下开采稳定分析中积累有益的经验。
1工程概况及采矿方法
1.1工程概况
五道羊岔矿床为太古宙岩浆晚期结晶分异形成的富铁矿床，矿区总面积5.006km2，开采标高为1060～-280m。槽探、钻探工程揭露显示，矿区共有10个矿组47条矿体，矿体剖面形态各异，呈隐伏、半隐伏型。
1.2开采范围及方法
五道羊岔矿区采区范围及相应的工程地质勘探线布置如图1。矿区分I～V共5个采区生产，共布置37条工程地质勘探线以查明地下矿体剖面形态。根据矿山现状及矿体赋存条件，确定了各采区的开采对象及开采方案如下。
1)I区：开采对象为11#～13#勘探线、17#～31#勘探线之间的矿体，开采标高1060～820m。自上而下分中段开采，采用分段空场法开采，水平方向采至矿体边界。
2)II区：开采对象为3#～8#勘探线、12#～32#勘探线之间的矿体，开采标高720～-280m。自上而下分中段开采，中段内采用自远而近后退式开采。采用房柱法开采，水平方向采至矿体边界。
3)III区：开采对象为45#～47#勘探线、51#～63#勘探线之间的矿体，开采标高1040～759m。自上而下分中段开采，中段内采用自远而近后退式开采。采用空场法开采，水平方向采至矿体边界。
4)IV区：开采对象为108#～102#勘探线之间的矿体，开采标高850～940m。整体自上而下分中段开采，中段内后退式开采。采用留矿法开采，水平方向采至矿体边界。
5)V区：开采对象为100#～110#勘探线之间的矿体，开采标高850～700m。整体自上而下分中段开采，中段内后退式开采。采用空场法开采，水平方向采至矿体边界。
2模拟计算
2.1计算模型
本次研究采用了FLAC2D二维有限差分数值模拟软件，对该矿区开采后的地表移动进行模拟。采用的计算为二维计算，计算平面为上述37条勘探线剖面。每条剖面的计算范围：水平方向上取沿该勘探线走向两端各延长1倍的距离，以消除边界效应的影响。对于每条计算剖面，采用摩尔-库伦模型模拟开挖计算。
模型两端水平方向上的位移固定，底部为固定边界，顶部为自由边界。模型上部自由边界上的网格节点为监测点。首先计算在没有开采状态下，由于自重应力引起的y向(垂直)沉降位移及速度，并使之归0。然后按照上述开采方案，模拟计算各剖面线上矿体开采完成后地表发生的沉降值及变形破坏范围。
2.2计算参数
计算剖面中涉及的岩性为铁矿体、斜长角闪岩、花岗岩。通过现场对不同岩性岩石的取样及室内岩石力学试验，结合相关工程经验类比进行经验折减，确定了计算力学参数(表1)。
2.3计算结果与分析
2.3.1地表动态沉陷特征
选取上述37条勘探线剖面，建立相应的数值计算模型。对每个剖面进行了仿真开采，分析开采完成后各条剖面上的位移场变化规律。由于篇幅所限，本节仅给出101#勘探线剖面仿真模拟开采的结果(图2)，旨在说明分析方法，其它可以此类推。不同开采阶段计算剖面上垂直位移等值线图如图2(a)和图2(c)。
可以看出在垂直方向上，随着开采的持续进行，地表沉降呈渐进式发展。不同开采阶段的地表沉陷位移曲线如图3。由图3可以看出，竖直位移曲线相对平滑，且连续性良好。在开采全部完成时地表沉降量达到最大，为9.7mm，出现在距勘探线左端400m的边坡处。不同开采阶段计算剖面的水平位移等值线图如图2(b)和图2(d)。
可见随着开采的进行，岩体水平位移的变化规律趋于明显，主要表现在采空区上下盘岩体向空区的挤近。即上、下盘岩体发生相向移动，上盘岩体总体上发生向左的移动，下盘岩体则发生向右的移动。沿勘探线走向的地表沉陷位移图如图3。
由图3可以看出，水平位移曲线近似以沉降中心为原点呈正弦曲线分布，最终的地表水平位移值为7.2mm，出现在距勘探线左端457m的边坡处。
2.3.2地表沉陷预测评价
通过对上述5个采区37条剖面不同开采阶段下的二维仿真计算，分别得出了各个剖面的垂直位移值和水平位移值。在上述计算剖面的离散化过程中，每条剖面计算模型上部自由边界上的网格节点为地表监测点。在开挖计算完成后，通过提取这些监测点的平面坐标及相应的位移量值，即可得到各剖面线上开采完成后相应监测点的空间三维数据。
基于Surfer平台，对上述的监测点数据进行处理。按《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》，以沉陷值为10mm作为沉陷盆地边界，绘制各采区的沉陷等值线图。整个矿区全部开采完成后的地表沉陷区划图如图4。
可以看出，地表沉降区主要分布在矿体上盘地表及下盘地表的一定区域内。I区、III区、IV区和V区地表呈椭圆形盆地塌陷形态，相应的各区最大沉降量分别为70、18、40、20mm。由于II区矿体厚度较大，在纵向延展较远，地表受开采影响较大，故开采结束后的地表沉降边界扩展较大，且无明显的沉降中心，最大沉降量值为40mm。
利用Surfer三维建模功能，对上述5矿区37个计算模型开采沉陷前后的二维面积在空间上的插值，可得到各矿区开采前后地表沉陷的体积变化，开采诱发的地表沉陷体积见表2。可以看出，由于II采区矿体规模大，在5个采区中II采区开采前后发生的塌陷量最大，约占整个矿区塌陷体积量的97.5%，这与上述数值模拟得到的结果是一致的。因此，应合理地布置II采区的矿山生产，以减少该采区矿山生产期间的运输成本及潜在矿山沉陷灾害带来的损失。
3结论
1)地表沉陷随着开采的进行呈渐进式发展，竖直位移曲线和水平位移曲线相对平滑;在对称性良好的水平开采情况下，竖直位移曲线的对称轴位置大致为正弦形态的水平位移分量的对称中心。
2)通过对不同采区地表位移量的分析和界定，给出了5个采区全部开采结束后的地表沉陷边界和地表塌陷体积。根据预测地表沉陷预测的位移量值，认为虽然II矿区地表建筑物在地表沉陷边界外，但其稳定性在一定程度上仍受开采沉陷的影响，在后续的矿山生产过程中应注意对这些建筑物的监测。
参考文献：
[1]陈元非，王磊，查剑锋.矿山开采沉陷预计与分析[J].煤矿安全，2017，48(1)：99-106.
[2]王军保，刘新荣，刘小军.开采沉陷动态预测模型[J].煤炭学报，2015，40(3)：516-521.
[3]崔希民，邓喀中.煤矿开采沉陷预计理论与方法研究评述[J].煤炭科学技术，2017，45(1)：160-169.
[4]刘宝琛，戴华阳.概率积分法的又来与研究进展[J].煤矿开采，2016，21(2)：1-3.
[5]朱晓峻，郭广礼，方齐.概率积分法预计参数反演方法研究进展[J].金属矿山，2015，44(4)：173-177.
[6]秦世界，张和生，李国栋.基于FLAC3D的煤矿开采沉陷预计及与概率积分法的对比分析[J].煤炭工程，2014，46(6)：96-102.
[7]范玮刚，冯光明，韩晓东，等.浅埋薄煤层超高水充填开采地表移动规律[J].能源技术与管理，2010(5)：19.
[8]刘剑，郝兵元，黄辉.大采宽工作面开采沉陷规律的FLAC3D模拟研究[J].煤炭技术，2014，33(11)：132.
[9]刘康，轩大洋.长壁全部充填与条带充填开采沉陷控制效果对比[J].煤矿安全，2013，44(11)：60-63.
[10]张文泉，刘海林，赵凯.厚松散层薄基岩条带开采地表沉陷影响因素研究[J].采矿与安全工程学报，2016，33(6)：1065-1071.
[11]刘玉成.煤层开采地表移动过程的FLAC3D模拟研究[J].煤炭科学技术，2012，40(5)：93-95.
推荐期刊：《煤炭科学技术》创刊于1973年，是由国家煤矿安全监察局主管、煤炭科学研究总院主办的综合性煤炭科技期刊。主要刊载煤炭工业各领域取得的最新成果、新技术、研究动态、新经验及新产品等方面的论文，内容包括煤矿开采、矿井建设、地质勘探、煤矿电气、煤矿机械、煤化工、科学管理、国外技术等。
