一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法

摘要

本发明公开了一种煤岩体分区注水渗流‑损伤‑应力耦合数值模拟方法，依据地质勘测结果，建立煤层模型；通过分区，分别对非采动影响区和采动影响区的煤岩体注水渗流进行模拟；在非采动影响区，编程通过计算比较各网格的拉剪力矩值，判断煤体是否在水压应力作用下发生形变甚至断裂；同时采用无网格法模拟煤体断裂过程，用边界元法模拟渗流过程，在微观尺度上结合两种模拟方法的优点；在采动影响区，分别采用N‑S方程和Darcy定律两种数学模型进行模拟，准确的模拟煤岩体分区注水过程中煤体损伤及水分的运移规律。为煤层注水提供坚实的数据基础，在大量减少了煤体被破碎为尘粒的可能性，降低了煤尘的产生量，保证了煤层开采的安全性。

说明书

技术领域

本发明涉及矿山岩石力学领域，尤其涉及一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法。

背景技术

煤层注水是通过钻孔，将压力水和水溶液注入煤体的过程，水的湿润作用使煤体塑性增强，脆性减弱，当煤体受外力作用时，许多脆性破碎变为塑性形变，因而大量减少了煤体被破碎为尘粒的可能性，降低了煤尘的产生量。整个过程涉及计算流体力学、断裂岩石力学与流固耦合等多门学科，对其进行数值模拟的过程较为复杂。

目前针对煤层注水的数值模拟研究大都以宏观角度进行，通过模拟宏观层面的渗流速度场及渗流压力场等研究水分在煤体内的运移规律。但作为一种典型的多孔介质材料，煤体润湿的本质是水分进入煤体内的众多细微孔隙，宏观角度的模拟无法复现这一过程以供深入研究。且现有技术中针对水分渗流的研究大都基于有限元或离散元之一单独进行分析，但是有限元难以实现煤体断裂，离散元无法准确描述渗流过程与水分增量数据。

因此，现有技术有待于更进一步的改进和发展。

发明内容

鉴于上述现有技术的不足，本发明的目的在于提供一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，更准确的模拟出煤岩体分区注水过程中煤体损伤及水分的运移规律，为煤层注水提供数据。

为解决上述技术问题，本发明方案包括：

一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其包括以下步骤：

步骤一，依据地质勘测结果，建立包含断层的煤岩体模型；

步骤二，对地质勘测得到之煤岩体中的相应煤体进行扫描，并结合FDK三维重建算法构建具有真实孔隙结构特征的三维数字岩心模型；

步骤三、在边界元环境下，执行“采动裂隙”生成算法，判断是否为采动影响区，即判断距离采掘面是否小于80m；如果是采动影响区，则执行步骤四；如果不是采动影响区，则执行步骤五；

步骤四、判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则进行N-S计算处理，储存并输出结果；若小于30μm²，则进行Darcy计算处理，储存并输出结果；

步骤五、判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则执行步骤六，若小于30μm²，则执行步骤十二；

步骤六、以C#语言为基础，编程以埋深、地质条件、上覆岩层岩性及煤层倾角生成初始应力集中点；

步骤七、对经过步骤六处理的三维数字岩心模型进行剖分网格；

步骤八、通过计算煤体网格所受拉剪力矩，将网格所受拉剪力矩与抗拉剪力矩比对，寻找最易被破坏的网格集合；若存在拉剪力矩大于抗拉剪力矩的网格，则标记为“原生裂隙”，然后执行步骤九；若不存在拉剪力矩大于抗拉剪力矩的网格，则终止计算，执行步骤十一；

步骤九、替换“原生裂隙”网格中的材料为煤层内气体；

步骤十、判断经过步骤九处理后之“原生裂隙”的总面积是否超过断层面积的两倍，如果没有超过，则生成删除“原生裂隙”网格后的新模型，返回重复步骤八；如果超过，则终止计算，执行步骤十一；

步骤十一、构建最终包含地质构造原生裂隙、计算过程中生成裂隙两类煤层裂隙的模型并重新剖分网格，导出STL通用几何后进行N-S计算处理，储存并输出结果；

步骤十二、在边界元环境下，耦合自定义方程计算渗流场与应力参数后，执行步骤十三；

步骤十三、在无网格仿真环境中，判断是否存在拉剪力矩大于煤岩体相关参数的点，如果有，则将此类点依次连接，闭合区域标记为“无效煤岩体”，执行步骤十四；如果没有，则执行步骤十五；

步骤十四、将步骤十三中，标记的“无效煤岩体”边界重新添加水分入口条件；

步骤十五、单独进行湍流模拟，并计算应力分布，并储存该节点的模拟结果；

步骤十六、判断步骤十五中模拟结果的断裂脱离是否达到煤体表面，如果没有，则执行步骤十七；如果到达煤体表面，则执行步骤十八；

步骤十七、累计存储时间是否达到预设模拟时间，如果没有达到，则返回执行步骤十二；如果达到计算时间，则执行步骤十八；

步骤十八、停止无网格仿真环境中的运算，仅在边界元环境中计算湍流与应力后则进行Darcy计算处理，储存并输出结果；

步骤十九、将步骤四、步骤十一与步骤十八的输出结果整合后输出并存储为独立文件，得到量化统计结果。

所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其中，N-S计算处理包括：

步骤A、N-S初始化；

步骤B、计算水压和瓦斯压；

步骤C、判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行步骤D，若小于，则直接执行步骤D；

步骤D、判断是否达到计算时间设定，若不能达到，则返回重新执行步骤A，若达到计算时间，则直接储存并输出结果。

所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其中，Darcy计算处理包括：

步骤E、Darcy初始化；

步骤F、计算水压和瓦斯压；

步骤J、判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行步骤H，若小于，则直接执行步骤H；

步骤H、判断是否达到计算时间设定，若不能达到，则返回重新执行步骤E，若达到计算时间，则直接储存并输出结果。

所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其中，上述步骤一具体的还包括：建立包含断层的煤岩体模型的过程中预留出“断层起点”、“断层终点”、“断层转折点”、“方位差”、“落差”与“倾角”六个参数供用户输入。

所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其中，上述步骤二聚体的还包括：利用图像处理算法对三维数字岩心模型进行滤波操作，平滑模型边缘，进而基于阈值实现数据分割得到煤体的微观孔隙结构，剔除连通性较差的孤岛孔隙，导出得到最终的STL格式煤体孔隙数字模型。

所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其中，上述量化统计结果包括水压场模拟结果、瓦斯压力模拟结果、煤体压力模拟结果、水分流速模拟结果与瓦斯流速模拟结果。

本发明提供了一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，分别采用Navier-Stokes方程和Darcy定律两种数学模型，基于边界元和无网格法两种方法进行模拟；边界元法与有限元相比，具有单元个数少，数据准备简单等优点，用边界元法模拟渗流过程，采用无网格法模拟煤体断裂过程，在微观尺度上结合两种模拟方法的优点，可以更准确的模拟出煤岩体分区注水过程中煤体损伤及水分的运移规律，为煤层注水提供坚实的数据基础，在大量减少了煤体被破碎为尘粒的可能性，降低了煤尘的产生量，保证了煤层开采的安全性。

附图说明

图1为本发明中煤岩体分区注水渗流-损伤-应力耦合数值模拟方法的示意图；

图2为本发明中N-S计算处理的示意图；

图3为本发明中Darcy计算处理的示意图；

图4为本发明中去除孤岛煤块示意图；

图5为本发明中去除孤岛后的煤体模型示意图；

图6为本发明中为水压结果模拟图；

图7为本发明瓦斯压力结果模拟图；

图8为本发明煤体压力结果模拟图；

图9为本发明水分流速结果模拟图；

图10为本发明瓦斯流速结果模拟图。

具体实施方式

本发明提供了一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提供了一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，如图1所示的，其包括以下步骤：

步骤一，依据地质勘测结果，建立包含断层的煤岩体模型；

步骤二，对地质勘测得到之煤岩体中的相应煤体进行扫描，并结合FDK三维重建算法构建具有真实孔隙结构特征的三维数字岩心模型；

步骤三、在边界元环境下，执行“采动裂隙”生成算法，判断是否为采动影响区，即判断距离采掘面是否小于80m；如果是采动影响区，则执行步骤四；如果不是采动影响区，则执行步骤五；

步骤四、判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，如图2所示的，则进行N-S计算处理，储存并输出结果；若小于30μm²，如图3所示的，则进行Darcy计算处理，储存并输出结果；

步骤五、判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则执行步骤六，若小于30μm²，则执行步骤十二；

步骤六、以C#语言为基础，编程以埋深、地质条件、上覆岩层岩性及煤层倾角生成初始应力集中点。

步骤七、对经过步骤六处理的三维数字岩心模型进行剖分网格；

步骤八、通过计算煤体网格所受拉剪力矩，将网格所受拉剪力矩与抗拉剪力矩比对，寻找最易被破坏的网格集合；若存在拉剪力矩大于抗拉剪力矩的网格，则标记为“原生裂隙”，然后执行步骤九；若不存在拉剪力矩大于抗拉剪力矩的网格，则终止计算，执行步骤十一；

步骤九、替换“原生裂隙”网格中的材料为煤层内气体；

步骤十、判断经过步骤九处理后之“原生裂隙”的总面积是否超过断层面积的两倍，如果没有超过，则生成删除“原生裂隙”网格后的新模型，返回重复步骤八；如果超过，则终止计算，执行步骤十一；

步骤十一、构建最终包含地质构造原生裂隙、计算过程中生成裂隙两类煤层裂隙的模型,并重新剖分网格，导出STL通用几何后进行N-S计算处理，储存并输出结果；

步骤十二、在边界元环境下，耦合自定义方程计算渗流场与应力参数后，执行步骤十三；

步骤十三、在无网格仿真环境中，判断是否存在拉剪力矩大于煤岩体相关参数的点，如果有，则将此类点依次连接，闭合区域标记为“无效煤岩体”，执行步骤十四；如果没有，则执行步骤十五；

步骤十四、将步骤十三中，标记的“无效煤岩体”边界重新添加水分入口条件；

步骤十五、单独进行湍流模拟，并计算应力分布，并储存该节点的模拟结果；

步骤十六、判断步骤十五中模拟结果的断裂脱离是否达到煤体表面，如果没有，则执行步骤十七；如果到达煤体表面，则执行步骤十八；

步骤十七、累计存储时间是否达到预设模拟时间，如果没有达到，则返回执行步骤十二；如果达到计算时间，则执行步骤十八；

步骤十八、停止无网格仿真环境中的运算，仅在边界元环境中计算湍流与应力后则进行Darcy计算处理，储存并输出结果；

步骤十九、将步骤四、步骤十一与步骤十八的输出结果整合后输出并存储为独立文件，得到量化统计结果。

在本发明的另一较佳实施例中，N-S计算处理包括：

步骤A、N-S初始化；

步骤B、计算水压和瓦斯压；

步骤C、判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行步骤D，若小于，则直接执行步骤D；

步骤D、判断是否达到计算时间设定，若不能达到，则返回重新执行步骤A，若达到计算时间，则直接储存并输出结果。

还有Darcy计算处理包括：

步骤E、Darcy初始化；

步骤F、计算水压和瓦斯压；

步骤J、判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行步骤H，若小于，则直接执行步骤H；

步骤H、判断是否达到计算时间设定，若不能达到，则返回重新执行步骤E，若达到计算时间，则直接储存并输出结果。

更进一步的，上述步骤一具体的还包括：建立包含断层的煤岩体模型的过程中预留出“断层起点”、“断层终点”、“断层转折点”、“方位差”、“落差”与“倾角”六个参数供用户输入。

并且上述步骤二聚体的还包括：利用图像处理算法对三维数字岩心模型进行滤波操作，平滑模型边缘，进而基于阈值实现数据分割得到煤体的微观孔隙结构，剔除连通性较差的孤岛孔隙，导出得到最终的STL格式煤体孔隙数字模型。

上述量化统计结果包括水压场模拟结果、瓦斯压力模拟结果、煤体压力模拟结果、水分流速模拟结果与瓦斯流速模拟结果。

为了更进一步的描述本发明，以下列举更为详尽的实施例。

一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法，其包括以下步骤：

步骤一、依据地质勘测结果，建立包含断层的煤岩体模型；

依据实际地质勘测的结果，编程建立煤层模型，其中，预留出“断层起点”、“断层终点”、“断层转折点”、“方位差”、“落差”、“倾角”六个参数供用户输入；

步骤二、CT阵列模型；

对地质勘测得到之煤岩体中的相应煤体进行扫描实验，结合FDK三维重建算法构建具有真实孔隙结构特征的三维数字岩心模型。利用图像处理算法对三维数字岩心模型进行滤波操作，平滑模型边缘，进而基于阈值实现数据分割得到煤体的微观孔隙结构，剔除连通性较差的孤岛孔隙，如图4中非浅色区域，处理后的煤体孔隙结构如图5所示，导出即为最终的STL格式煤体孔隙数字模型。

步骤三、在边界元环境下，执行“采动裂隙”生成算法，判断是否为采动影响区，即判断距离采掘面是否小于80m。如果是采动影响区，则执行步骤四；如果不是采动影响区，则执行步骤五。

步骤四、判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则执行步骤六，若小于30μm²，则执行步骤七。

步骤五、判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则执行步骤八，若小于30μm²，则执行步骤十四。

步骤六、进行N-S计算，储存并输出结果。

步骤七、进行Darcy计算，储存并输出结果。

步骤八、以C#语言为基础，编程以埋深、地质条件、上覆岩层岩性及煤层倾角生成初始应力集中点。

步骤九、剖分网格。

步骤十、通过计算煤体网格所受拉剪力矩，将网格所受拉剪力矩与抗拉剪力矩比对，寻找最易被破坏的网格集合。若存在拉剪力矩大于抗拉剪力矩的网格，则标记为“原生裂隙”，然后执行步骤十一。若不存在拉剪力矩大于抗拉剪力矩的网格，则终止计算，执行步骤十三。

步骤十一、替换“原生裂隙”网格中的材料为煤层内气体。

步骤十二、判断“原生裂隙”的总面积是否超过断层面积的两倍，如果没有超过，则生成删除原生裂隙网格的新模型，返回重复步骤十；如果超过，则终止计算，执行步骤十三。

步骤十三、构建最终包含地质构造原生裂隙、计算过程中生成裂隙两类煤层裂隙的模型，并重新剖分网格，导出STL通用几何后执行步骤六。

步骤十四、在边界元环境下，耦合自定义方程计算渗流场与应力参数后，执行步骤十五。

步骤十五、在无网格仿真环境中，编写相应算法，判断是否存在拉剪力矩大于煤岩体相关参数的点，如果有，则将此类点连接，闭合区域标记为“无效煤岩体”，执行步骤十六；如果没有，则执行步骤十七。

步骤十六、将步骤十五中，标记的“无效煤岩体”边界重新添加水分入口条件。

步骤十七、单独进行湍流模拟，并计算应力分布，并储存该节点的模拟结果。

步骤十八、判断断裂脱离是否达到煤体表面，如果没有，则执行步骤十九；如果到达煤体表面，则执行步骤二十。

步骤十九、累计存储时间是否达到预设模拟时间，如果没有达到，则返回执行步骤十四；如果达到计算时间，则执行步骤二十。

步骤二十、停止无网格法环境中的运算，仅在边界元环境中计算湍流与应力后执行步骤七。

步骤二十一、将整合后的结果输出并存储为独立文件，量化统计结果，如图6的水压场模拟结果图、图7的瓦斯压力模拟结果图、图8的煤体压力结果模拟图、图9的水分流速模拟结果图与图10的瓦斯流速模拟结果图。

还包括N-S计算处理和Darcy计算处理的步骤：

N-S计算处理步骤；

步骤六一、N-S初始化。

步骤六二、计算水压和瓦斯压。

步骤六三、判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行步骤六四，若小于，则直接执行步骤六四。

步骤六四、判断是否达到计算时间设定，若不能达到，则返回重新执行步骤六一，若达到计算时间，则直接储存并输出结果。

Darcy计算处理步骤；

步骤七一、Darcy初始化。

步骤七二、计算水压和瓦斯压。

步骤七三、判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行步骤七四，若小于，则直接执行步骤七四。

步骤七四、判断是否达到计算时间设定，若不能达到，则返回重新执行步骤七一，若达到计算时间，则直接储存并输出结果。

实施例1

Step 1依据地质勘测结果，建立包含断层的煤岩体模型，其具体包括如下步骤：

(1)在笛卡尔坐标系下，依据注水煤层高度、走向长度及倾向长度建立煤体几何模型，模型维度为三维；

(2)其中，预留出“断层起点”、“断层终点”、“断层转折点”、“方位差”、“落差”、“倾角”六个参数供用户输入。

(3)然后导出STL通用几何。

Step 2CT阵列模型，其具体包括如下步骤：

(1)采用精度0.5μm的nanoVoxel-2000系列-X射线三维显微镜对长焰煤进行扫描实验，结合FDK三维重建算法构建了具有真实孔隙结构特征的三维数字岩心模型；

(2)利用图像处理算法对三维数字岩心模型进行滤波操作，减少噪声，平滑模型边缘，然后基于灰度值，进行体渲染，进而基于阈值实现数据分割得到煤体的微观孔隙结构；

(3)剔除连通性较差的孤岛孔隙，如图4中非浅蓝色区域(具体参见其他证明材料中的图4)，处理后的煤体孔隙结构如图5所示(具体参见其他证明材料中的图5)，导出即为最终的STL格式煤体孔隙数字模型

Step 3在边界元环境下，执行“采动裂隙”生成算法，判断是否为采动影响区，即距离采掘面是否小于80m。如果是采动影响区，则执行步骤四；如果不是采动影响区，则执行步骤五。为了理论本身的完备性，引入以下假设：

(1)忽略或不考虑基质岩体的透水性和储水性(因为相比较于“采动裂隙”，基质岩体的透水性和储水性能均特别弱)。

(2)“采动裂隙”渗流服从达西定律其中k＝-b²/12

(3)“采动裂隙”变形规律服从Goodman节理模型。

Step 4判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则执行步骤六，若小于30μm²，则执行步骤七。

Step 5判断流域的截面积是否大于30μm²，若流域截面积大于30μm²，则执行步骤八，若小于30μm²，则执行步骤十四。

Step 6进行N-S计算，储存并输出结果。

(1)N-S初始化。

以N-S来描述水分注水孔及裂隙中水分的自由流动，在笛卡尔坐标系下表示为：

式中，t—时间；ρ—煤体中的水压，MPa；v—运动黏滞系数；u_x、u_y、u_z—分别为x、y、z轴轴向单位质量力，mg·s^-2；F_x、F_y、F_z是外力的分量，N；▽—那勃勒算子；u、v、w是流体在t时刻，在点(x、y、z)处的速度分量。

(2)计算水压和瓦斯压。

(3)判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行(4)，若小于，则直接执行(4)。

(4)判断是否达到计算时间设定，若不能达到，则返回重新执行(1)，若达到计算时间，则直接储存并输出结果。

Step 7进行Darcy计算，储存并输出结果。

(1)Darcy初始化。

(2)计算水压和瓦斯压。

Darcy定律描述水在煤体内的渗流运动，其微分方程在笛卡尔坐标系中表示为：

式中，P为煤体中的水压MPa，V_x与V_y分别为x、y轴方向上的速度分量；K为渗透率，μ为水的动力粘度Pa·s,g为重力加速度常数；ρ为液体密度。

Darcy定律认为煤层内煤层气运动基本符合线性渗透规律，即：

其中，V为煤层气运动的速度矢量m/s；gradp为煤层内煤层气孔隙压力梯度Pa/m；μ为煤层内煤层气孔隙压力函数；k为煤层气渗透性系数为煤层气的绝对黏度,Pa·s；

瓦斯沿煤层的流动一般服从Darcy定律，即瓦斯的渗流速度与瓦斯压力梯度成正比：

其中，q为瓦斯渗流速度，cm/s；λ为煤层透气系数，m²/(MPa²·d)；μ为流体的动力黏度系数,Pa·s；k为渗透率，m²；为瓦斯压力梯度，P/cm。

(3)判断水压是否大于瓦斯压，若大于，则将下一网格材料变更为水，执行(4)，若小于，则直接执行(4)。

(4)判断是否达到计算时间设定，若不能达到，则返回重新执行(1)，若达到计算时间，则直接储存并输出结果。

Step 8以C#语言为基础，编程以埋深、地质条件、上覆岩层岩性及煤层倾角生成初始应力集中点。

Step 9剖分网格，其具体包括如下步骤：

对煤体几何模型进行四面体网格剖分，在计算机运算力允许的范围内尽可能提高网格剖分精度，以保证裂隙计算结果的准确性与可靠性。

Step 10边界元环境下，计算并比较各网格的拉剪力矩值，其具体包括如下步骤：

(1)依据但不限于Coulomb-Mohr强度准则描述岩体的拉剪破裂过程。Coulomb准则认为岩体的拉剪破裂过程发生在某一平面产生的破坏拉剪力矩值超过了抗拉剪力矩值即材料的内聚力和乘以常数的平面法应力。其数学表达式为|τ|＝C+γtanμ其中τ是拉剪力矩值；C是内聚力或者黏结力，是无正压力时的抗拉剪强度；μ是内摩擦角；γ是一个固定常数。

(2)通过计算煤体网格所受拉剪力矩，将网格所受拉剪力矩与抗拉剪力矩比对，寻找最易被破坏的网格集合。若存在拉剪力矩大于抗拉剪力矩的网格，则标记为“采动裂隙”，然后执行Step 11。若不存在拉剪力矩大于抗拉剪力矩的网格，则终止计算，执行Step 13。

Step 11替换“原生裂隙”网格中的材料为煤层内气体。

(1)在煤体几何模型中剔除Step 10中被标记为“采动裂隙”的网格；

(2)将被剔除的几何单元材料(网格)变更为煤层内气体，孔隙率变更为1；

(3)对因该单元被剔除而新生成的边界设定与相邻边界相同的边界条件；

Step 12判断“原生裂隙”的总面积是否超过断层面积的两倍，其具体包括如下步骤：

(1)如果没有超过，则生成删除原生裂隙网格的新模型，返回重复Step 10；(2)如果超过，则终止计算，执行Step 13。

Step 13构建最终包含地质构造原生裂隙、计算过程中生成裂隙两类煤层裂隙的模型,其具体包括如下步骤:

(1)生成包含断层的煤岩体模型；

(2)重新剖分网格；

(3)导出STL通用几何；

Step 14在边界元环境下，耦合自定义方程计算渗流场与应力参数。

(1)依据但不限于达西定律描述水分在煤体中的运移过程，即以作为渗流模拟的数学模型，其中，t为时间；φk为裂隙比空隙度；s为沿裂隙长度的坐标；k_f为裂隙渗透系数；W为源汇项；p为煤体中水压(MPa)。若研究人员有需要添加的其它流动控制方程，可手动添加；

(2)岩体在孔隙流体的作用下，遵循修正的应力规律。依据但不限于基于流体渗流压力的应力方程σ_ij＝σ'_ij+αpσ_ij计算煤体应力分布，其中α称为有效应力系数。实践证明，α是孔隙压p和体积应力θ的函数，α＝f(p,θ)。

在包含其它应力分布规律或特殊情况的煤层，可编程加入其它本构方程，对更加贴近实际的煤体应力分布进行模拟分析。

Step 15在无网格仿真环境中，编写算法，判断是否存在拉剪力矩大于煤岩体相关参数的点，其具体包括如下步骤:

(1)如果存在，将拉剪力矩大于煤岩体相关参数的点，连接起来，形成闭合区域标记为“无效煤岩体”，继续执行Step 16；

(2)如果没有拉剪力矩大于煤岩体相关参数的点，则继续执行Step 17。

Step 16将Step 15中，标记的“无效煤岩体”边界重新添加水分入口条件。

Step 17单独进行湍流模拟，并计算应力分布，并储存该节点的模拟结果。

Step 18判断断裂脱离是否达到煤体表面，其具体包括如下步骤:

(1)如果断裂脱离没有到达煤体表面，则执行Step 19；

(2)如果断裂脱离到达煤体表面，则执行Step 20。

Step 19累计存储时间是否达到预设模拟时间，其具体包括如下步骤:

(1)如果存储时间没有达到预设模拟时间，则返回执行Step 14；

(2)如果存储时间达到预设模拟时间，则执行Step 20。

Step 20停止无网格法环境中的运算，仅在边界元环境中计算湍流与应力，后执行Step 7。

Step 21将整合后的结果输出并存储为独立文件，量化统计结果，如图6的水压场模拟结果图、图7的瓦斯压力模拟结果图、图8的煤体压力结果模拟图、图9的水分流速模拟结果图与图10的瓦斯流速模拟结果图。

依本数值模拟方法结果运算得到的结果可以看出，随着更大范围的煤体被水分润湿，水分运移的主要动力来源已变为毛细作用力而非注水压力，因此水分与瓦斯的平均压力均逐渐下降；但水分压力的下降梯度较为恒定，而瓦斯压力收到不规则孔、裂隙结构的影响，下降梯度差异明显。煤体随注水作业进行而逐渐被软化，因此综合作用压力总体上呈现出下降趋势。同时，水分与瓦斯的流速变化与各自的压力变化基本一致，但瓦斯流速高于水分流速。

上述现象与工程现场应用结果相一致，说明本发明的模拟方法所得结果是可靠性的。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。