一种考虑垮落带采空区的计算机数值模拟方法与流程

本发明属于煤炭采出后采空区受上覆岩层压缩响应的计算机数值模拟技术领域，具体是一种考虑垮落带采空区的计算机数值模拟方法。

背景技术：

煤炭从地下采出引起上覆岩层移动和应力重新分布。长壁开采为我国主要的采煤方法，开采后形成采空区，采空区是采场煤炭采出的产物，对地形地貌影响巨大，采空区响应对覆岩运移影响显著。

岩层移动的关键层理论，通过试验、理论研究和现场证据证明了长壁开采覆岩形成三带，即垮落带、裂隙带和弯曲下沉带。图1示为长壁采场覆岩运动规律。垮落带岩体呈散体和碎裂结构，岩块间空隙大、连通性好；裂隙带岩层断裂失去连续性，层状特性仍存在；弯曲下沉带只弯曲不断裂，保持连续性和原始特征。为垮落带采空区发育特征提供了基础依据。

许多地面工程需要确定采空区响应，除了地上工程，许多地下工程位于采空区或其附近，同样涉及采空区响应，如分层开采上分层采空区下布置下分层巷道、沿空留巷、沿空掘巷以及错层位开采技术等。

采空区除了影响地形地貌外，还对采动应力分布具有决定性的重要作用，鉴于采空区的重要性，国内学者对采空区材料及其响应进行了研究。

大量对采空区破碎矸石压实特性的实验研究表明，采空区矸石压缩应力-应变曲线呈指数或近似指数的增长趋势。其中，作为该类试验的权威代表，pappas等通过在不同煤矿不同长壁工作面取矸石样品进行多组实验室压缩试验，对目前提出的多个矸石压缩应力-应变公式分析，得出salamon总结出的矸石压缩的应力-应变曲线最符合矸石的压实特性。

然而，采矿科学在煤矿设计中常需要进行相关数值模拟以研究设计的合理性和科学性。如，数值模拟技术对巷道应力的分布、煤柱的设计如今已应用非常广泛和流行。科学的数值模拟方法将大幅提升模拟的准确性，从而提升模拟结论的现实指导性。尽管目前垮落带采空区的实验室试验已较多，获得较多的数据和结论，但目前流行的长壁开采数值模拟均对采空区不做任何处理，及在模型中挖空后便进行运算，采空区的支撑作用没有任何体现，即覆岩压力没有传递至采空区，也没有通过采空区传递至底板，模拟结果必然与实际大相径庭。垮落带采空区的数值模拟技术的，没有一种科学的模拟方法，限制了计算机技术在采矿工程领域的指导作用。目前计算机的运算速度和容量以迅雷不及掩耳之势发展，计算机更新和换代令人应接不暇，利用好计算机，提出科学合理的采空区数值模拟方法必能大幅提高计算机在采矿领域的应用价值和指导意义。

技术实现要素：

本发明为了解决上述问题，提供一种考虑垮落带采空区的计算机数值模拟方法。

本发明采取以下技术方案：一种考虑垮落带采空区的计算机数值模拟方法，所述方法包括以下步骤：

第一步：flac3d建模；

第二步：确定除采空区以外的岩体力学参数；

第三步：将采空区材料赋予双屈服本构模型；

第四步：确定采空区材料相关参数；

第五步：将所有岩体力学和采空区参数代入flac3d模型中进行运算，便可得到考虑采空区的数值模拟结果。

在所述的第一步中，建模尺寸数据由钻孔地质资料和工作面尺寸确定，其中，flac3d模型中煤层岩层的高度及岩性与钻孔获得的岩芯相对应，模型中工作面尺寸按照实际中工作面的尺寸划定；通过地应力测量、模型埋深以及质地构造，确定模型的边界条件，边界条件包括模型四周的固定情况、测压系数(水平应力与竖直应力的比值)以及顶部加载情况。处理方式为：顶部边界自由，底面限制竖直方向的运动，前后左右四个面限制水平运动，顶部可自由运动；顶部加载根据上方未建模型的高度确定，若未建高度为h，覆岩容重为γ，则顶部施加竖直向下的力为γh；侧压系数则根据实测水平应力与竖直应力的比值在模型中赋值。

在所述的第二步中，对岩样进行实验室岩石力学试验，获得rocdata软件需要的实验室数据来估算σci,mi.和a并得到广义h-b强度参数mb,s和a，并得到最适m-c准则等价摩擦角和粘聚力c，这些参数用于模型中各岩层的赋值和运算，选取除采空区以外的岩体材料的本构模型，确定各岩层的物理力学参数，本构模型表达式为：

式中，σ1为破坏时的最大主应力；σ3为破坏时的最小主应力或是三轴试验中的围压；mb为岩体的h-b常数；s为与岩体特性有关的材料常数，反映岩体破碎程度，其取值范围0～1；a为表征节理岩体的常数；σci为完整岩石的单轴抗压强度ucs；而：

式中，d为节理岩体弱化因子，主要考虑爆破破坏和应力松弛对节理岩体的扰动程度，取值为0～1；gsi为地质强度指标，gsi的值由工程岩体的岩体结构、结构面特征因素综合确定；

h-b与mohr-coulomb准则，下文简写为h-c准则，拟合确定最适m-c准则等价摩擦角和粘聚力c的表达式为：

式中，σ3n＝σ3max/σci，σ3max为最小主应力上限值，

抗压强度和抗拉强度分别为:

σc＝σcis^a

要得到岩体的广义h-b强度参数mb,s和a,需要各岩层的以下参数:完整岩石的单轴抗压强度σci；完整岩石参数mi；地质强度指标gsi；扰动系数d，通过实验室数据拟合或通过软件内置典型岩层参数范围表估算σci,mi.和a并得到广义h-b强度参数mb,s和a以及最适m-c准则等价摩擦角和粘聚力c，其中σci、mi、gsi和d参数可通过实验室数据拟合，缺少实验室数据的情况下可根据软件内置典型岩石经验参数得到σci,mi和a，该内置参数通过世界范围内的大量岩石工程收集所得。

在所述的第三步中，采空区材料赋予双屈服本构模型，该本构面模型中应力-应变呈类似指数式增长的方式。该本构模型适于模拟压应力引起体积永久缩减的低胶结粒状散体岩土材料。模型中考虑了由于各向同性压力引起的永久性体积变化，还考虑了flac程序中应变软化/硬化模型中剪切破坏和拉破坏判断准则以及体积屈服面。

在所述的第四步中，双屈服模型中所需的参数有：体积模量、剪切模量、剪胀角和摩擦角参数。

参数通过在flac3d中进行采空区材料压缩试验与salamon经验公式迭代匹配得到；

采空区矸石应力-应变曲线符合salamon经验公式：

式中，σ为加载于材料上的轴向压力，[mpa]；ε为在该压力下的应变，m/m；e0为初始切线模量[mpa]；εm为最大可能应变；且：

其中，σc为垮落区矸石单轴抗压强度，εm只与初始碎胀系数b有关：

初始碎胀系数为：

其中，hc为煤层上方垮落带高度；m为煤厚。

在所述的第四步中，在flac3d中建立简单块体模型，底部固支，顶部不约束，其它四个面约束横向位移，通过在其表面施加竖向速率模拟加载过程，通过迭代试错法来确定体积模量、剪切模量、剪胀角和摩擦角参数，经过反复试验，可最终确定参数，使参数满足数值模拟所得的应力-应变曲线与salamon公式曲线的匹配效果理想。

与现有技术相比，本发明提供了一种考虑垮落带采空区的计算机数值模拟方法。该方法为了改善目前数值模拟中不考虑采空区造成模拟结果不符合实际情况的现状。不考虑采空区，采空区不发挥任何支撑作用，即覆岩压力没有传递至采空区，也没有通过采空区传递至底板，则会造成支承压力偏大，覆岩和底板破环和运移范围偏大，本发明与现广泛应用的技术相比具有以下优点：模型单元随机生成，而非常规矩形，避免人为设定网格造成的影响；煤岩体参数则基于h-b强度准则和m-c准则拟合得出，考虑了岩体结构的弱面，模拟数据更加可靠；将垮落角考虑入数值模拟，更加符合实际。采用双屈服模型模拟采空区,其参数基于迭代试错匹配与salamon应力-应变曲线拟合确定，实验室试验数据与计算机模拟数据的匹配确保模拟结果的可靠性。

附图说明

图1为覆岩三带发育示意图；

图2为flac3d数值模型示例；

图3为rocdata软件h-b准则与m-c准则换算实例；

图4为双屈服本构模型表征的力学响应图像；

图5为相似模拟所得垮落角示意图；

图6为数值模拟与salamon经验公式曲线匹配结果。

具体实施方式

一种考虑垮落带采空区的计算机数值模拟方法，包括以下步骤：

第一步：建模。建模对计算的影响显著，建议采用外接建模接口软件，如ansys等进行建模，可建立出flac3d自身无法建立的随机分布的网格模型，计算更加可靠。根据钻孔数据、工作面尺寸等逐层建立数值模型。通过地应力测试，确定模型的边界条件。图2所示为一flac3d模型。

第二步：选取除采空区以外的岩体材料的本构模型，确定各岩层的物理力学参数。准确的本构模型是得出可靠模拟结果的基础。hoek-brown经验强度准则(下文简称h-b准则)是由hoeke和brownet，在大量矿山岩体力学测试、室内实验基础并根据大量工程实践经验得出，该本构模型考虑了岩体的不连续性，已在世界范围内得到广泛应用和实践检验。该本构模型表达式为：

式中，σ1为破坏时的最大主应力；σ3为破坏时的最小主应力(或是三轴试验中的围压)；mb为岩体的h-b常数；s为与岩体特性有关的材料常数，反映岩体破碎程度，其取值范围0～1；a为表征节理岩体的常数；σci为完整岩石的单轴抗压强度(ucs)。而：

式中，d为节理岩体弱化因子，主要考虑爆破破坏和应力松弛对节理岩体的扰动程度，取值为0～1；gsi为地质强度指标，gsi的值由工程岩体的岩体结构、结构面特征等因素综合确定。

h-b与mohr-coulomb准则(下文简写为h-c准则)拟合确定最适m-c准则等价摩擦角和粘聚力c的表达式为：

式中，σ3n＝σ3max/σci，σ3max为最小主应力上限值，与特定的工程类型有关。

抗压强度和抗拉强度分别为:

σc＝σcis^a

要得到岩体的广义h-b强度参数(mb,s和a),需要各岩层的以下参数:完整岩石的单轴抗压强度σci；完整岩石参数mi；地质强度指标gsi；扰动系数d。rocdata能够通过实验室数据拟合或通过软件内置典型岩层参数范围表估算σci,mi.和a并得到广义h-b强度参数mb,s和a以及最适m-c准则等价摩擦角和粘聚力c。其中σci、mi、gsi和d参数可通过实验室数据拟合，缺少实验室数据的情况下可根据软件内置典型岩石经验参数得到σci,mi和a，该内置参数通过世界范围内的大量岩石工程收集所得。这些参数将用于模型中各岩层的赋值和运算。图3所示为一换算过程例子。

第三步：确定采空区材料的本构模型。合理处理采空区才能得到合理的模拟结论。采空区矸石在逐渐压实过程中强度和模量逐渐增大，双屈服本构模型适用于模拟压应力引起体积永久缩减的低胶结粒状散体岩土材料，符合采空区材料的力学行为。除了将剪切和拉伸破坏包络线纳入flac外，通过将体积屈服面即应变软化/硬化模型(或称为“帽盖”)纳入flac3d，双屈服模型将各向同性的压力引起的永久体积改变考虑在内。与采空区具有类似性质的材料在压缩时除了发生剪切屈服外，还会发生一种和剪切屈服明显相互作用的屈服模式，例如液压回填材料或者由混凝土胶结形成的小颗粒材料，而双屈服本构模型能够表征这种行为。采空区材料赋予双屈服本构模型的目的是使采空区材料满足在逐渐压实过程中强度和模量呈salamon经验公式的指数式增长，符合试验现象，与实际接近。图4所示为双屈服破坏模型的力学响应。

第四步：确定采空区材料相关参数。双屈服模型中所需的参数有：体积模量、剪切模量、剪胀角和摩擦角参数。这些参数则通过在flac3d中进行采空区材料压缩试验与salamon经验公式迭代匹配得到。

采空区矸石应力-应变曲线符合salamon经验公式：

式中，σ为加载于材料上的轴向压力，[mpa]；ε为在该压力下的应变，m/m；e0为初始切线模量[mpa]；εm为最大可能应变。且：

其中，σc为垮落区矸石单轴抗压强度。εm只与初始碎胀系数b有关：

初始碎胀系数为：

其中，hc为煤层上方垮落带高度，；m为煤厚。

由此可见，确定采空材料参数，需要确定垮落带采空区范围。有若干种确定采空区范围的方式，包括经验公式、相似模拟实验、穿层钻孔、地面钻孔或其它地球物理勘探手段，以及上述两种或几中方法的结合。图5所示为相似模拟所得垮落角示意图。注意，相似模拟得出的垮落带范围的准确性不如现场实测得出的可靠，视工程精度要求适当选择采空区范围测定方法。

确定采空区范围(包括垮落带高度、长度、宽度以及垮落角)后，便可通过以上公式计算出salamon经验公式的参数e0和εm，salamon经验公式便可确定。

在flac3d中建立简单块体模型，底部固支，顶部不约束，其它四个面约束横向位移。通过在其表面施加竖向速率模拟加载过程。通过迭代试错法来确定体积模量、剪切模量、剪胀角和摩擦角参数。经过反复试验，可最终确定参数，使参数满足数值模拟所得的应力-应变曲线与salamon公式曲线的匹配效果理想。图6所示为匹配结果。

第五步：将所有岩体力学和采空区参数代入模型进行运算，得到考虑采空区的数值模拟结果。