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

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

Info

Publication number
CN109063257B
CN109063257B CN201810705645.3A CN201810705645A CN109063257B CN 109063257 B CN109063257 B CN 109063257B CN 201810705645 A CN201810705645 A CN 201810705645A CN 109063257 B CN109063257 B CN 109063257B
Authority
CN
China
Prior art keywords
coal
result
rock
grid
stress
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810705645.3A
Other languages
English (en)
Other versions
CN109063257A (zh
Inventor
周刚
尹文婧
张文政
魏星
张国宝
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong University of Science and Technology
Original Assignee
Shandong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong University of Science and Technology filed Critical Shandong University of Science and Technology
Priority to CN201810705645.3A priority Critical patent/CN109063257B/zh
Priority to RU2020110457A priority patent/RU2743121C1/ru
Priority to PCT/CN2018/100462 priority patent/WO2020006818A1/zh
Publication of CN109063257A publication Critical patent/CN109063257A/zh
Application granted granted Critical
Publication of CN109063257B publication Critical patent/CN109063257B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • EFIXED CONSTRUCTIONS
    • E21EARTH OR ROCK DRILLING; MINING
    • E21BEARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B49/00Testing the nature of borehole walls; Formation testing; Methods or apparatus for obtaining samples of soil or well fluids, specially adapted to earth drilling or wells
    • E21B49/006Measuring wall stresses in the borehole
    • EFIXED CONSTRUCTIONS
    • E21EARTH OR ROCK DRILLING; MINING
    • E21CMINING OR QUARRYING
    • E21C41/00Methods of underground or surface mining; Layouts therefor
    • E21C41/16Methods of underground mining; Layouts therefor
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mining & Mineral Resources (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geochemistry & Mineralogy (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Environmental & Geological Engineering (AREA)
  • Fluid Mechanics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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

Description

一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法
技术领域
本发明涉及矿山岩石力学领域,尤其涉及一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法。
背景技术
煤层注水是通过钻孔,将压力水和水溶液注入煤体的过程,水的湿润作用使煤体塑性增强,脆性减弱,当煤体受外力作用时,许多脆性破碎变为塑性形变,因而大量减少了煤体被破碎为尘粒的可能性,降低了煤尘的产生量。整个过程涉及计算流体力学、断裂岩石力学与流固耦合等多门学科,对其进行数值模拟的过程较为复杂。
目前针对煤层注水的数值模拟研究大都以宏观角度进行,通过模拟宏观层面的渗流速度场及渗流压力场等研究水分在煤体内的运移规律。但作为一种典型的多孔介质材料,煤体润湿的本质是水分进入煤体内的众多细微孔隙,宏观角度的模拟无法复现这一过程以供深入研究。且现有技术中针对水分渗流的研究大都基于有限元或离散元之一单独进行分析,但是有限元难以实现煤体断裂,离散元无法准确描述渗流过程与水分增量数据。
因此,现有技术有待于更进一步的改进和发展。
发明内容
鉴于上述现有技术的不足,本发明的目的在于提供一种含地质构造的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,更准确的模拟出煤岩体分区注水过程中煤体损伤及水分的运移规律,为煤层注水提供数据。
为解决上述技术问题,本发明方案包括:
一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其包括以下步骤:
步骤一,依据地质勘测结果,建立包含断层的煤岩体模型;
步骤二,对地质勘测得到之煤岩体中的相应煤体进行扫描,并结合FDK三维重建算法构建具有真实孔隙结构特征的三维数字岩心模型;
步骤三、在边界元环境下,执行“采动裂隙”生成算法,判断是否为采动影响区,即判断距离采掘面是否小于80m;如果是采动影响区,则执行步骤四;如果不是采动影响区,则执行步骤五;
步骤四、判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则进行N-S计算处理,储存并输出结果;若小于30μm2,则进行Darcy计算处理,储存并输出结果;
步骤五、判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则执行步骤六,若小于30μm2,则执行步骤十二;
步骤六、以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μm2,若流域截面积大于30μm2,如图2所示的,则进行N-S计算处理,储存并输出结果;若小于30μm2,如图3所示的,则进行Darcy计算处理,储存并输出结果;
步骤五、判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则执行步骤六,若小于30μm2,则执行步骤十二;
步骤六、以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μm2,若流域截面积大于30μm2,则执行步骤六,若小于30μm2,则执行步骤七。
步骤五、判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则执行步骤八,若小于30μm2,则执行步骤十四。
步骤六、进行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=-b2/12
(3)“采动裂隙”变形规律服从Goodman节理模型。
Step 4判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则执行步骤六,若小于30μm2,则执行步骤七。
Step 5判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则执行步骤八,若小于30μm2,则执行步骤十四。
Step 6进行N-S计算,储存并输出结果。
(1)N-S初始化。
以N-S来描述水分注水孔及裂隙中水分的自由流动,在笛卡尔坐标系下表示为:
式中,t—时间;ρ—煤体中的水压,MPa;v—运动黏滞系数;ux、uy、uz—分别为x、y、z轴轴向单位质量力,mg·s-2;Fx、Fy、Fz是外力的分量,N;—那勃勒算子;u、v、w是流体在t时刻,在点(x、y、z)处的速度分量。
(2)计算水压和瓦斯压。
(3)判断水压是否大于瓦斯压,若大于,则将下一网格材料变更为水,执行(4),若小于,则直接执行(4)。
(4)判断是否达到计算时间设定,若不能达到,则返回重新执行(1),若达到计算时间,则直接储存并输出结果。
Step 7进行Darcy计算,储存并输出结果。
(1)Darcy初始化。
(2)计算水压和瓦斯压。
Darcy定律描述水在煤体内的渗流运动,其微分方程在笛卡尔坐标系中表示为:
式中,P为煤体中的水压MPa,Vx与Vy分别为x、y轴方向上的速度分量;K为渗透率,μ为水的动力粘度Pa·s,g为重力加速度常数;ρ为液体密度。
Darcy定律认为煤层内煤层气运动基本符合线性渗透规律,即:
其中,V为煤层气运动的速度矢量m/s;gradp为煤层内煤层气孔隙压力梯度Pa/m;μ为煤层内煤层气孔隙压力函数;k为煤层气渗透性系数为煤层气的绝对黏度,Pa·s;
瓦斯沿煤层的流动一般服从Darcy定律,即瓦斯的渗流速度与瓦斯压力梯度成正比:
其中,q为瓦斯渗流速度,cm/s;λ为煤层透气系数,m2/(MPa2·d);μ为流体的动力黏度系数,Pa·s;k为渗透率,m2为瓦斯压力梯度,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为沿裂隙长度的坐标;kf为裂隙渗透系数;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的瓦斯流速模拟结果图。
依本数值模拟方法结果运算得到的结果可以看出,随着更大范围的煤体被水分润湿,水分运移的主要动力来源已变为毛细作用力而非注水压力,因此水分与瓦斯的平均压力均逐渐下降;但水分压力的下降梯度较为恒定,而瓦斯压力收到不规则孔、裂隙结构的影响,下降梯度差异明显。煤体随注水作业进行而逐渐被软化,因此综合作用压力总体上呈现出下降趋势。同时,水分与瓦斯的流速变化与各自的压力变化基本一致,但瓦斯流速高于水分流速。
上述现象与工程现场应用结果相一致,说明本发明的模拟方法所得结果是可靠性的。
当然,以上说明仅仅为本发明的较佳实施例,本发明并不限于列举上述实施例,应当说明的是,任何熟悉本领域的技术人员在本说明书的教导下,所做出的所有等同替代、明显变形形式,均落在本说明书的实质范围之内,理应受到本发明的保护。

Claims (6)

1.一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其包括以下步骤:
步骤一,依据地质勘测结果,建立包含断层的煤岩体模型;
步骤二,对地质勘测得到之煤岩体中的相应煤体进行扫描,并结合FDK三维重建算法构建具有真实孔隙结构特征的三维数字岩心模型;
步骤三、在边界元环境下,执行“采动裂隙”生成算法,判断是否为采动影响区,即判断距离采掘面是否小于80m;如果是采动影响区,则执行步骤四;如果不是采动影响区,则执行步骤五;
步骤四、判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则进行N-S计算处理,储存并输出结果;若小于30μm2,则进行Darcy计算处理,储存并输出结果;
步骤五、判断流域的截面积是否大于30μm2,若流域截面积大于30μm2,则执行步骤六,若小于30μm2,则执行步骤十二;
步骤六、以C#语言为基础,编程以埋深、地质条件、上覆岩层岩性及煤层倾角生成初始应力集中点;
步骤七、对经过步骤六处理的三维数字岩心模型进行剖分网格;
步骤八、通过计算煤体网格所受拉剪力矩,将网格所受拉剪力矩与抗拉剪力矩比对,寻找最易被破坏的网格集合;若存在拉剪力矩大于抗拉剪力矩的网格,则标记为“原生裂隙”,然后执行步骤九;若不存在拉剪力矩大于抗拉剪力矩的网格,则终止计算,执行步骤十一;
步骤九、替换“原生裂隙”网格中的材料为煤层内气体;
步骤十、判断经过步骤九处理后的“原生裂隙”的总面积是否超过断层面积的两倍,如果没有超过,则生成删除“原生裂隙”网格后的新模型,返回重复步骤八;如果超过,则终止计算,执行步骤十一;
步骤十一、构建最终包含地质构造原生裂隙、计算过程中生成裂隙两类煤层裂隙的模型,并重新剖分网格,导出STL通用几何后进行N-S计算处理,储存并输出结果;
步骤十二、在边界元环境下,耦合自定义方程计算渗流场与应力参数后,执行步骤十三;
步骤十三、在无网格仿真环境中,判断是否存在拉剪力矩大于煤岩体相关参数的点,如果有,则将此类点依次连接,闭合区域标记为“无效煤岩体”,执行步骤十四;如果没有,则执行步骤十五;
步骤十四、将步骤十三中,标记的“无效煤岩体”边界重新添加水分入口条件;
步骤十五、单独进行湍流模拟,并计算应力分布,并储存节点的模拟结果;
步骤十六、判断步骤十五中模拟结果的断裂脱离是否达到煤体表面,如果没有,则执行步骤十七;如果到达煤体表面,则执行步骤十八;
步骤十七、累计存储时间是否达到预设模拟时间,如果没有达到,则返回执行步骤十二;如果达到计算时间,则执行步骤十八;
步骤十八、停止无网格仿真环境中的运算,仅在边界元环境中计算湍流与应力后则进行Darcy计算处理,储存并输出结果;
步骤十九、将步骤四、步骤十一与步骤十八的输出结果整合后输出并存储为独立文件,得到量化统计结果。
2.根据权利要求1所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其特征在于,N-S计算处理包括:
步骤A、N-S初始化;
步骤B、计算水压和瓦斯压;
步骤C、判断水压是否大于瓦斯压,若大于,则将下一网格材料变更为水,执行步骤D,若小于,则直接执行步骤D;
步骤D、判断是否达到计算时间设定,若不能达到,则返回重新执行步骤A,若达到计算时间,则直接储存并输出结果。
3.根据权利要求1所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其特征在于,Darcy计算处理包括:
步骤E、Darcy初始化;
步骤F、计算水压和瓦斯压;
步骤J、判断水压是否大于瓦斯压,若大于,则将下一网格材料变更为水,执行步骤H,若小于,则直接执行步骤H;
步骤H、判断是否达到计算时间设定,若不能达到,则返回重新执行步骤E,若达到计算时间,则直接储存并输出结果。
4.根据权利要求1所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其特征在于,上述步骤一具体的还包括:建立包含断层的煤岩体模型的过程中预留出“断层起点”、“断层终点”、“断层转折点”、“方位差”、“落差”与“倾角”六个参数供用户输入。
5.根据权利要求1所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其特征在于,上述步骤二聚体的还包括:利用图像处理算法对三维数字岩心模型进行滤波操作,平滑模型边缘,进而基于阈值实现数据分割得到煤体的微观孔隙结构,剔除连通性较差的孤岛孔隙,导出得到最终的STL格式煤体孔隙数字模型。
6.根据权利要求1所述的煤岩体分区注水渗流-损伤-应力耦合数值模拟方法,其特征在于,上述量化统计结果包括水压场模拟结果、瓦斯压力模拟结果、煤体压力模拟结果、水分流速模拟结果与瓦斯流速模拟结果。
CN201810705645.3A 2018-07-02 2018-07-02 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法 Active CN109063257B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201810705645.3A CN109063257B (zh) 2018-07-02 2018-07-02 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法
RU2020110457A RU2743121C1 (ru) 2018-07-02 2018-08-14 Способ численного моделирования связи фильтрации/повреждений/напряжений при впрыске воды в каменноугольный массив при районировании
PCT/CN2018/100462 WO2020006818A1 (zh) 2018-07-02 2018-08-14 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810705645.3A CN109063257B (zh) 2018-07-02 2018-07-02 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法

Publications (2)

Publication Number Publication Date
CN109063257A CN109063257A (zh) 2018-12-21
CN109063257B true CN109063257B (zh) 2019-04-26

Family

ID=64818148

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810705645.3A Active CN109063257B (zh) 2018-07-02 2018-07-02 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法

Country Status (3)

Country Link
CN (1) CN109063257B (zh)
RU (1) RU2743121C1 (zh)
WO (1) WO2020006818A1 (zh)

Families Citing this family (39)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109117589B (zh) * 2018-09-11 2022-10-04 中煤科工集团重庆研究院有限公司 一种煤层顶板裂隙场的定量化描述方法
CN109883920A (zh) * 2019-03-08 2019-06-14 西南石油大学 表征岩石热损伤的方法及装置
CN110043259B (zh) * 2019-04-18 2020-05-26 中国矿业大学 一种大采深大采高复合顶工作面分区域回采工艺
CN110118991B (zh) * 2019-05-16 2020-06-23 中国矿业大学 一种基于微震损伤重构的采动应力评估方法
CN110489832B (zh) * 2019-07-31 2023-05-23 中国航发沈阳发动机研究所 一种用于湍流控制屏单元体气动性能的仿真试验方法
CN110390176B (zh) * 2019-07-31 2020-05-19 西南交通大学 一种无砟轨道冻结与损伤行为计算方法
CN110872942B (zh) * 2019-10-24 2021-09-28 中国石油化工股份有限公司 油藏注采耦合方式下不同作用力对采收率贡献的研究方法
CN110866337B (zh) * 2019-11-12 2021-06-01 中南大学 一种基于差应力的采动断层活化倾向性判别方法
CN111307687B (zh) * 2020-03-05 2022-07-19 中煤科工集团重庆研究院有限公司 高分子材料与煤岩体粘结的渗透性及强度评价方法
CN111666699B (zh) * 2020-04-30 2023-06-02 山东大学 基于rev全区域覆盖的岩体工程跨尺度模拟计算方法
CN111860952B (zh) * 2020-06-16 2024-02-27 重庆大学 一种突出煤层关键开采参数快速优选方法
CN111695285B (zh) * 2020-06-17 2023-12-22 大连海事大学 一种各向异性岩体应力-损伤-渗流耦合数值模拟方法
CN111898187B (zh) * 2020-07-23 2022-10-11 武汉大学 隧洞开挖渗流模拟分析的纵向模型范围取值的确定方法
CN112131709B (zh) * 2020-08-25 2024-04-19 山东大学 基于近场动力学本构模型的节理岩体力学仿真方法及系统
CN112507575A (zh) * 2020-11-03 2021-03-16 辽宁工程技术大学 基于flac 3d数值模拟的断层注浆加固效果评价方法
CN112329312B (zh) * 2020-11-10 2022-07-26 河海大学 一种三维渗流应力耦合内聚力单元的快速生成方法
CN112560226B (zh) * 2020-11-25 2023-12-26 国家能源投资集团有限责任公司 煤矿地下水库环境风险评估方法
CN112924331B (zh) * 2021-01-12 2022-10-04 江苏师范大学 水溶液浸泡后煤岩抗压强度的水岩耦合模型的建立方法
CN112800597B (zh) * 2021-01-13 2022-09-13 安徽马钢张庄矿业有限责任公司 一种矿山资源高中段智能精细高效生态开采分析方法
CN112903966A (zh) * 2021-01-20 2021-06-04 中国矿业大学(北京) 基于能量传递守恒的煤矿开采损伤范围确定方法
CN112926270B (zh) * 2021-03-18 2023-05-05 西安科技大学 瓦斯多因素耦合关系分析及预警模型构建方法
CN113533042B (zh) * 2021-07-07 2022-04-05 北京科技大学 一种表征岩石应力与破裂的综合性指标计算方法及应用
CN113654477B (zh) * 2021-08-16 2023-02-21 中国矿业大学 一种煤体变形测试装置、测试系统及测试方法
CN113722960A (zh) * 2021-08-31 2021-11-30 中国地质大学(武汉) 一种湿度运移和干缩开裂的三维数值模拟方法
CN113742940B (zh) * 2021-09-16 2023-09-08 重庆大学 一种采动卸压边界时空曲线簇确定方法
CN114002129B (zh) * 2021-12-02 2022-08-09 河北省交通规划设计研究院有限公司 一种高水压裂隙岩体渗流试验平台
CN114217048B (zh) * 2021-12-10 2024-05-24 国家能源投资集团有限责任公司 采矿三维模拟实验模型实验方法
CN114495679B (zh) * 2022-01-25 2022-10-28 中国矿业大学 一种真实煤二维微流体模型制作方法
CN114544437B (zh) * 2022-02-28 2023-12-19 中国矿业大学 一种煤岩流场原位荧光菌显微示踪方法
CN114660268B (zh) * 2022-03-22 2023-06-20 中铁水利水电规划设计集团有限公司 水库淹没区的抬田区耕地保水层渗流监测系统
CN114638137B (zh) * 2022-03-31 2023-03-28 王永亮 一种基于热-流-固-损伤耦合的干热岩产热预测方法
CN114809992B (zh) * 2022-04-20 2023-08-08 太原理工大学 一种低渗储层煤系气全生命周期高效抽采方法
CN115201898B (zh) * 2022-05-20 2023-06-20 中国地震局地球物理研究所 三维的注采诱发地震断层破裂滑动的数值模拟方法及系统
CN115290528B (zh) * 2022-07-04 2024-08-02 中国矿业大学 一种富水岩层裂隙结构-渗流特性相互作用的表征方法
CN115014982B (zh) * 2022-08-09 2022-10-11 中国矿业大学(北京) 一种基于气体运移压力波动特征评价煤岩损伤劣度的方法
CN115524275B (zh) * 2022-10-18 2023-08-22 中国水利水电科学研究院 考虑岩块渗流和裂隙渗流的裂隙岩体渗透张量确定方法
CN116756985B (zh) * 2022-11-29 2024-01-30 华东师范大学 基于COMSOL Multiphysics的场地多介质环境有机污染物运移模拟方法
CN117408191B (zh) * 2023-12-15 2024-04-02 山东大学 一种裂隙岩体渗透-流动注浆模拟方法及系统
CN117627669A (zh) * 2024-01-26 2024-03-01 中交一航局第三工程有限公司 一种基于盲区超前导洞扩挖的矩形顶管施工方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SU796420A1 (ru) * 1977-04-04 1981-01-15 Всесоюзный Научно-Исследовательскийинститут Горной Геомеханики Имаркшейдерского Дела Способ имитации напр женного сос-ТО Ни гОРНОгО МАССиВА HA МОдЕл Х изэКВиВАлЕНТНыХ МАТЕРиАлОВ
US8898046B2 (en) * 2010-12-10 2014-11-25 Baker Hughes Incorporated Method to improve reservoir simulation and recovery from fractured reservoirs
CN102252957B (zh) * 2011-04-19 2013-07-03 河南理工大学 卸压构造煤固-流转化参数实验测定装置及方法
CN103034765B (zh) * 2012-12-14 2015-08-05 天津大学 基于数值模拟的采空区注浆动态全过程仿真方法
CN104653226B (zh) * 2014-12-26 2017-07-14 中国矿业大学 一种基于应力梯度的煤矿冲击地压危险区域的划分方法
CN105401939A (zh) * 2015-11-30 2016-03-16 中国石油大学(北京) 一种多因素耦合作用下的煤层井壁稳定性分析方法
CN105787220B (zh) * 2016-04-22 2018-10-12 山东科技大学 一种煤层高压注水致裂-渗流数值模拟方法
CN106960070B (zh) * 2016-12-28 2020-02-21 山东科技大学 一种基于有限元-离散元ct重构煤体的渗流模拟方法
CN108038282A (zh) * 2017-11-30 2018-05-15 安徽理工大学 一种松散承压含水层下开采载荷传递数值模拟方法

Also Published As

Publication number Publication date
CN109063257A (zh) 2018-12-21
WO2020006818A1 (zh) 2020-01-09
RU2743121C1 (ru) 2021-02-15

Similar Documents

Publication Publication Date Title
CN109063257B (zh) 一种煤岩体分区注水渗流-损伤-应力耦合数值模拟方法
Liu et al. Study on the effect of cemented natural fractures on hydraulic fracture propagation in volcanic reservoirs
CN106960070B (zh) 一种基于有限元-离散元ct重构煤体的渗流模拟方法
Pine et al. The development of a new numerical modelling approach for naturally fractured rock masses
WO2016192077A1 (zh) 一种致密气压裂水平井数值试井模型建立求解方法
CN101446196B (zh) 三重介质油藏分支水平井的试井分析方法及装置
CN107060746A (zh) 一种复杂裂缝性油藏流动模拟的方法
CN106991244B (zh) 一种基于图论的裂隙网络连通性及渗流计算的方法
CN113901681A (zh) 一种全寿命周期页岩气储层双甜点三维可压性评估方法
Lin et al. A finite-discrete element based appoach for modelling the hydraulic fracturing of rocks with irregular inclusions
CN106886046B (zh) 确定缝洞型气藏未投产区块可动用储量的方法
CN107423466A (zh) 一种支撑剂嵌入和裂缝导流能力定量预测的数值模拟方法
CN108386176A (zh) 一种天然裂缝与人工裂缝延伸规律物模试验方法
CN106649963A (zh) 体积压裂复杂缝网平均裂缝长度和等效裂缝条数确定方法
CN114462272B (zh) 一种深层复杂构造下页岩气水平井井眼轨迹优化方法
Chen et al. Effect of dominated coal pores and fractures on water migration after low-pressure water injection based on CT images
LI et al. Research on random propagation method of hydraulic fracture based on zero-thickness cohesive element
Yang et al. Influence of reservoirs/interlayers thickness on hydraulic fracture propagation laws in low-permeability layered rocks
CN118261075A (zh) 一种水平井重复压裂增能暂堵压裂参数设计方法
CN112100890B (zh) 一种压裂水平缝起裂扩展全三维数学模型的计算方法
CN105929461A (zh) 一种动静态岩石力学参数校正系统
Butlanska Cone penetration test in a virtual calibration chamber
CN115705454A (zh) 一种基于相场法的裂缝扩展模拟压裂设计优化方法
CN108643894A (zh) 三维油藏物理模型断层设置方法
Yang et al. Influence of the different fluids on the hydraulic fracturing mechanism of sandy mudstone

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant