特厚煤层综放开采工作面底板破坏深度的连续-离散耦合模拟方法

技术领域

本发明涉及底板突水机理与防治、煤炭绿色开采技术领域，具体涉及一种特厚煤层综放开采工作面底板破坏深度的连续-离散耦合模拟方法。

背景技术

晋陕蒙宁甘赋存我国近60％的煤炭资源，但水资源占有量仅为全国的8％。更严重的是，地下煤炭开采会对该地区珍贵的奥陶系灰岩含水层造成严重破坏，导致人畜饮水困难，严重制约当地经济社会发展。煤系地层由不同岩性的完整岩块、大量裂隙构成。底板内丰富的天然裂隙，以及采动引起的裂隙，是导致底板突水与水资源破坏的决定性因素之一。煤矿生产实践表明：特厚煤层综放开采工作面底板的破坏深度往往大于中厚煤层。采用合理的数值模拟方法预测采动底板，对于优化特厚煤层综放工作面参数、绿色采煤至关重要。

近年来，国内外学者主要采用有限元、拉格朗日元等数值程序开展采动底板破坏深度研究。有限元、拉格朗日元均基于连续介质假设而建立，采用弹性、塑性、损伤力学理论，研究采动过程中底板的应力、变形、塑性区深度与演化规律。但是，采动底板岩石往往出现破裂，而不是塑性变形；底板中大量发育的天然裂隙随机分布、且尺度不一，在数值建模中需要生成随机分布的天然裂隙；更为重要的是，特厚煤层综放开采条件下，底板岩石先后承担准静态压缩-准静态卸载-动力压缩，只有有效模拟出来这种特殊荷载，方能数值预测采动底板裂隙发育深度。

发明内容

本发明的目的在于，提供一种特厚煤层综放开采工作面底板破坏深度的连续-离散耦合模拟方法。

为了实现上述任务，本发明采取如下的技术解决方案：

一种特厚煤层综放开采工作面底板破坏深度的连续-离散耦合模拟方法，其特征在于，包括以下步骤：

S1：确定底板天然裂隙间距，生成包含天然裂隙的几何模型；

S2：对几何模型划分单元，并在天然裂隙位置处嵌入裂隙单元；

S3：实验确定岩石压缩力学参数，以及I型、Ⅱ型断裂力学参数，岩石剪切摩擦参数；

S4：编制VUSDFLD子程序，自动删除煤集合单元，模拟工作面推进过程。

S5：施加数值模型边界条件，输入材料力学参数，设置增量步大小，设置输出参量，包括，应力、应变、位移、单元节点力、裂隙单元的状态变量。

S6：通过连续-离散耦合模拟程序，求解模型，得到特厚煤层综放工作面底板裂隙的时空演化规律。

S7：采动底板导水通道演化分析。根据数值模拟结果，对超前支承压力、采空区底板卸压、特厚煤层综放动压这种循环加卸载条件下的底板导水通道发育规律进行分析。进而结合地下含水层分布情况，优化煤炭开采方案，实现绿色采煤。

根据本发明，在步骤S1中，采用钻孔窥视方法和钻孔取芯所得各岩层岩芯的平均长度，确定煤层及围岩内主要天然裂隙的方位、沿钻孔轴向方向的裂隙间距。在此实测数据基础上，通过Matlab生成一组满足x,y,z三个正交方向上满足裂隙间距的空间离散点。而后利用Delaunay三角形生成相应的Voronoi多面体。将多面体导入Part模块，生成与实测间距、方位结果相近的天然裂隙几何模型。

在步骤S2中，包括如下内容：

S2.1、对S1几何模型进行网格划分，生成实体单元；

S2.2、在实体单元边界处嵌入裂缝单元。

其中，步骤S2.1中，对S1几何模型进行网格划分，生成实体单元的内容是：

在Mesh模块下，采用Seed Part模式，对S1中建立的几何体进行布种操作。布置种子的间距不超过Voronoi多面体棱长的1/2。通过Assign Mesh Controls控制单元的形状为六面体、四面体或者以六面体为主。而后通过Mesh Part对整个模型进行网格划分。进一步地，在Assign Element Type按钮下，将划分好的实体单元修改为C3D8实体单元类型。

步骤S2.2中，在实体单元边界处嵌入裂缝单元的内容是：

(1)对C3D8实体单元的节点编号进行更新。在完成步骤S2.1之后，在Job模块下创建一个新的作业，并在相应工作目录下找到对应作业的.inp文件。打开.inp文件，在关键字*nset和*elset下方数字中，找到最大节点和最大单元编号，并将其分别记为nset

对于C3D8实体单元，任一单元(假设编号为elset

在关键字“*Part”之后“*Assembly”之前添加关键字*Elset＝solid，而后按照修改后的单元编号顺序写入相应数字。从而得到实体单元集合。

(2)创建裂隙单元节点和单元编号。在.inp文件中，复制单元编号不同但有2个或者4个节点编号相同的数组。对于这种数组，记相邻两个单元的单元编号为e

而裂缝单元的编号根据实体单元最大编号elset

将初始.inp文件另存为，并在新的.inp文件中添加关键字*Element type＝coh3d8，同时在关键字“*Part”之后“*Assembly”之前添加关键字*Elset＝coh，而后按照修改后的单元、节点编号顺序写入相应数字。从而在实体单元边界处嵌入裂隙单元，并得到裂隙单元集合。

在步骤S3中，包括如下内容：

S3.1、实验确定岩石压缩力学参数；

S3.2、实验确定I型、Ⅱ型断裂力学参数；

S3.3、实验确定岩石剪切摩擦参数。

步骤S3.1中，实验确定岩石压缩力学参数的内容是：

试验室温度控制在(21±0.5)℃，以减小温度对试验数据采集设备的扰动。将直径50mm高度100mm的岩石试件放置在三轴压缩试验机中，用位移控制方式，加载速率为0.02mm/min，考虑了0MPa、5MPa和10MPa共3种围压条件。加载过程分为2个阶段，首先通过压力室内的液压油施加围压至预定围压值，然后通过轴向位移控制形成偏应力直至岩样破坏。记录压缩破坏过程中岩石应力-应变曲线，计算弹性模量、强度。

步骤S3.2中，实验确定I型、Ⅱ型断裂力学参数的内容是：

对煤岩试件开展紧凑拉伸实验，确定I性断裂参数。煤岩试件中最小几何尺寸至少为煤岩内晶体尺寸的十倍。建议实验所用煤岩试件为扁圆柱体，直径为100mm，高度为35mm。用0.5mm厚度的金刚石切割片沿着试件对称轴预制一条裂隙，其长度为15mm，两个对称加载孔直径为13mm；加载孔圆心到试件直边的距离为13.0mm，到试件对称中心的距离为8mm。通过加载辊轮将试件与试验机相连接。在实验过程中，位移加载速率为0.02mm/min，记录试件变形、破坏过程中荷载-裂纹尖端张开位移曲线，得到I型断裂强度、断裂能、断裂位移w

对煤岩试件开展剪切贯穿实验，确定Ⅱ性断裂参数。剪切贯穿试件为圆柱体，直径为50mm，高度为50mm。沿圆柱轴线方向在试件两平行断面相对预制两条不同深度的同心圆裂缝，对中心处的完整圆柱形区域施加载荷，预制裂缝半径25mm，宽度1.5mm，上部预制裂缝深度10mm，下部预制裂缝深度30mm。进行贯穿剪切实验时，采用直径为25mm的圆柱形加载钢头对煤岩试件中心完整小圆柱部分进行施加载荷，采用位移加载方式，加载速度为0.02mm/min。实验过程中记录荷载-位移曲线。由此计算得到Ⅱ型断裂强度、断裂能、断裂位移。

步骤S3.3中，实验确定岩石剪切摩擦参数的内容是：

制备圆柱体试件，直径50mm，高度100mm。而后用沿圆柱体轴向方向巴西劈裂，制备含单一粗糙裂隙的试件。将试件置于三轴压力室中，将围压提高至0.2MP，确保所有压力、位移传感器正常。根据地应力实测数据将围压、轴压同步增至试验设定值。而后采用位移加载方式，以恒定速率0.002mm/s进行剪切位移加载，设置剪切位移为实验设备最大加载位移的1/2。实验过程中记录剪切荷载-剪切位移曲线。由此计算得到剪切强度、剪切模量。

在步骤S4中，包括如下内容：

S4.1：在前处理程序中设置需要删除的单元集合；

S4.2：编制VUSDFLD子程序；

S4.3：将子程序导入主程序，模拟工作面推进过程；

其中，步骤S4.1中，在前处理程序中设置需要删除的单元集合的内容是：

在步骤S2的数值模型中设置需要删除的单元集合，具体如下：

(1)首先确定数值模型中需要删除煤层单元的范围；

(2)对步骤S2中建立的单元集合“*Elset＝solid”、“*Elset＝coh”分别复制，然后通过减选方式框定需要删除的煤层单元；

(3)根据实体单元平均尺寸，以及工作面实际推进速度，确定每一次推进距离下需要删除单元的数量。然后通过减选方式，进一步确定每一次推进距离下需要删除的实体单元和裂隙单元，得到各自的单元编号。

步骤S4.2中，编制VUSDFLD子程序的内容是：

第一，在编译器中，首先定义导入VUSDFLD子程序的参数。采用“subroutinevusdfld”子程序，导入主程序中的参数，包括单元编号、积分点坐标、状态变量参数、应力、位移、总时间、分析步时间、时间增量、方向等，并定义各参数为标量、矢量及张量。

第二，定义每一开挖步需要删除的单元集合。在定义好总时间totalTime、每推进一定距离所需要的时间deleteTime、单元号nblock、单元积分点坐标coordMp、方向direct的基础上，固定需要删除的起始单元积分点坐标coordMp(0)。在此基础上，当deleteTime与totalTime的比值为c时(c为某个常数)，定义所有满足该条件的单元为“集合1”；当deleteTime与totalTime的比值为2c时，定义所有满足该条件的单元为“集合2”；依此类推。

第三，设置初始时刻所有单元的初始状态变量stateOld＝1。而后通过stateNew＝0命令，依次删除“集合1”、“集合2”等单元集合。由此实现单元删除、工作面推进的目标。

步骤S4.3中，将子程序导入主程序的内容是：

在Job模块下，新建立一个作业任务，在“通用”按钮下选择子程序文件，即可将步骤S4.2的VUSDFLD.for子程序导入主程序。

在步骤S6中，连续-离散耦合模拟程序的内容是：

连续-离散耦合模拟程序思路为：在初始边界条件、和地应力条件下，通过有限元方法计算数值模型所有节点的节点力。通过步骤S4.2的VUSDFLD.for子程序删除单元后，采用有限元方法计算实体单元节点力，而后通过实体单元与裂隙单元的共享节点传递应力，利用有限元方法计算裂隙单元的位移量w。采用离散元方法计算裂隙面的剪切摩擦效应。实体单元(模拟完整煤岩块体)、裂隙单元(模拟非贯通裂隙)、接触对(模拟贯通裂隙)的力学参数根据步骤S3确定。通过这种连续-离散耦合模拟程序，求解工作面推进过程中底板岩体应力、应变、位移等规律。

采动底板裂隙生成的数学判定准则为：根据上述连续-离散耦合模拟程序，计算裂隙单元的节点力、节点位移。通过下式计算任一时刻裂隙单元的断裂能：

其中，G

当有G

本发明的特厚煤层综放开采工作面底板破坏深度的连续-离散耦合模拟方法，与传统的有限元、拉格朗日元等基于连续介质假设的数值模拟方法相比，本发明能够实时监测采动底板裂隙发育规律、重现采动底板承受超前支承压力加载-采空区卸载-覆岩运动所致动载，进而提供适用于特厚煤层综放开采工作面底板导水裂隙演化的模拟方法。带来的技术创新在于：

1、特厚煤层综放开采工作面，底板承受的是静力-动力的加/卸载条件，具体而言，是承受超前支承压力加载-采空区卸载-覆岩运动导致的动载。本发明的连续-离散耦合模拟方法，可有效模拟这种特殊载荷。是对传统有限元、离散元、拉格朗日元方法的重大改进。

2、为地下煤层开采引起底板破坏、裂隙时空演化规律研究提供了有力的数值模拟工具。可数值实现底板岩石在采动应力状态下从连续介质转化为离散介质的全过程，能够实时监测采动底板裂隙发育规律。

3、可用于分析不同开采参数、不同地质参数影响下采动底板导水裂隙时空演化，从而优化开采参数、实现地下煤炭绿色开采、为进一步底板注浆改造提供依据。

附图说明

图1是连续-离散耦合数值模拟流程示意图；

图2是数值计算模型及相应的裂隙单元集合；

图3是I型断裂试件；

图4是Ⅱ型断裂试件；

图5是采动底板裂隙分布；

图6是特厚煤层综放开采工作面开采条件下围岩位移场；

图7是特厚煤层综放开采工作面开采条件下围岩应力场。

下面结合附图和实施例对本发明做进一步地详细说明说明。

具体实施方式

本发明的设计思路是，针对煤层底板裂隙随机分布，采用基于钻孔窥视方法和钻孔取芯所得各岩层岩芯的平均长度数据，建立裂隙随机分布的几何模型；对几何模型划分实体单元，并在天然裂隙与实体单元边界位置处嵌入裂隙单元；通过实验获取岩石压缩力学参数，I型、Ⅱ型断裂力学参数，以及岩石剪切摩擦参数；编制反映采煤工作面推进过程的Vusdfld子程序，将其导入连续-离散耦合模拟主程序，数值实现特厚煤层综放开采工作面底板从完整、碎裂、结构面剪切滑移的全过程。

本实施例给出一种特厚煤层综放工作面矿压显现的有限离散元模拟方法，具体包括以下步骤：

S1：确定底板天然裂隙间距，生成包含天然裂隙的几何模型；

本实施例中，采用钻孔窥视方法和钻孔取芯所得各岩层岩芯的平均长度，确定煤层及围岩内主要天然裂隙的方位、沿钻孔轴向方向的裂隙间距；在此实测数据基础上，通过Matlab生成一组满足x,y,z三个正交方向上满足裂隙间距的空间离散点；而后利用Delaunay三角形生成相应的Voronoi多面体；将多面体导入Part模块，生成与实测间距、方位结果相近的天然裂隙几何模型。

通过钻孔窥视方法和钻孔取芯，确定各类型岩层中主要结构面的间距如图2所示。通过Voronoi多面体生成程序，建立长x＝200m，高y＝174m，厚度z＝1m的三维几何模型，如图2所示。

S2：对几何模型划分单元，并在天然裂隙位置处嵌入裂隙单元；

具体包括如下内容：

S2.1：对S1几何模型进行网格划分，生成实体单元，具体内容是：

在Mesh模块下，采用Seed Part模式，对S1中建立的几何体进行布种操作；布置种子的间距不超过Voronoi多面体棱长的1/2；通过Assign Mesh Controls控制单元的形状为六面体、四面体或者以六面体为主，而后通过Mesh Part对整个模型进行网格划分。进一步地，在Assign Element Type按钮下，将划分好的实体单元修改为C3D8实体单元类型。

S2.2：在实体单元边界处嵌入裂缝单元，具体内容是：

(1)对C3D8实体单元的节点编号进行更新

在完成步骤S2.1之后，在Job模块下创建一个新的作业，并在相应工作目录下找到对应作业的.inp文件；打开.inp文件，在关键字*nset和*elset下方数字中，找到最大节点和最大单元编号，并将其分别记为nset

对于C3D8实体单元，假设编号为elset

在关键字“*Part”之后“*Assembly”之前，添加关键字*Elset＝solid，而后按照修改后的单元编号顺序写入相应数字，从而得到实体单元集合。

(2)创建裂隙单元节点和单元编号

在.inp文件中，复制单元编号不同但有2个或者4个节点编号相同的数组；对于这种数组，记相邻两个单元的单元编号为e

而裂缝单元的编号根据实体单元最大编号elset

将初始.inp文件另存，并在新的.inp文件中添加关键字*Element type＝coh3d8，同时在关键字“*Part”之后“*Assembly”之前添加关键字*Elset＝coh，而后按照修改后的单元、节点编号顺序写入相应数字；从而在实体单元边界处嵌入裂隙单元，并得到裂隙单元集合。

在Mesh模块下，采用Seed Part模式，对几何模型进行布种操作。布置种子的间距为3m。通过Assign Mesh Controls控制单元的形状为六面体，采用Mesh Part对整个模型进行网格划分。而后通过Assign Element Type，将划分好的实体单元修改为C3D8实体单元类型。将所有实体单元组成名为“solid”的单元集合。

在实体单元边界处嵌入裂隙单元，单元类型为coh3d8，并将所有裂隙单元组合成名为“coh”的单元集合(图2中折线所示)。

S3：实验确定岩石压缩力学参数，以及I型、Ⅱ型断裂力学参数，岩石剪切摩擦参数；

具体包括如下内容：

S3.1：实验确定岩石压缩力学参数

试验室温度控制在21℃±0.5℃，以减小温度对试验数据采集设备的扰动，将直径50mm高度100mm的岩石试件放置在三轴压缩试验机中，用位移控制方式，加载速率为0.02mm/min，考虑0MPa、5MPa和10MPa共3种围压条件，加载过程分为2个阶段；首先通过压力室内的液压油施加围压至预定围压值，然后通过轴向位移控制形成偏应力直至岩样破坏，记录压缩破坏过程中岩石应力-应变曲线，计算弹性模量、强度；

S3.2：实验确定I型、Ⅱ型断裂力学参数

对煤岩试件开展紧凑拉伸实验，确定I性断裂参数：

煤岩试件中最小几何尺寸至少为煤岩内晶体尺寸的十倍；实验所用煤岩试件为扁圆柱体，直径为100mm，高度为35mm；用0.5mm厚度的金刚石切割片沿着试件对称轴预制一条裂隙，其长度为15mm，两个对称加载孔直径为13mm；加载孔圆心到试件直边的距离为13.0mm，到试件对称中心的距离为8mm；通过加载辊轮将试件与试验机相连接；

在实验过程中，位移加载速率为0.02mm/min，记录试件变形、破坏过程中荷载-裂纹尖端张开位移曲线，得到I型断裂强度、断裂能、断裂位移w

对煤岩试件开展剪切贯穿实验，确定Ⅱ性断裂参数：

剪切贯穿试件为圆柱体，直径为50mm，高度为50mm，沿圆柱轴线方向在试件两平行断面相对预制两条不同深度的同心圆裂缝，对中心处的完整圆柱形区域施加载荷，预制裂缝半径25mm，宽度1.5mm，上部预制裂缝深度10mm，下部预制裂缝深度30mm；进行贯穿剪切实验时，采用直径为25mm的圆柱形加载钢头对煤岩试件中心完整小圆柱部分进行施加载荷，采用位移加载方式，加载速度为0.02mm/min；

实验过程中记录荷载-位移曲线，由此计算得到Ⅱ型断裂强度、断裂能、断裂位移；

S3.3：实验确定岩石剪切摩擦参数

制备圆柱体试件，直径50mm，高度100mm。而后用沿圆柱体轴向方向巴西劈裂，制备含单一粗糙裂隙的试件；将试件置于三轴压力室中，将围压提高至0.2MP，确保所有压力、位移传感器正常；根据地应力实测数据将围压、轴压同步增至试验设定值，而后采用位移加载方式，以恒定速率0.002mm/s进行剪切位移加载，设置剪切位移为实验设备最大加载位移的1/2；

实验过程中记录剪切荷载-剪切位移曲线。由此计算得到剪切强度、剪切模量。

本实施例中，对底板典型岩层取芯，并制备成标准圆柱体试件(直径50mm，高度100mm)、圆柱体试件(直径50mm，高度200mm并预制裂隙)、扁圆柱体试件(直径50mm，高度50mm并预制环形裂隙)，如图3和图4所示。分别通过压缩试验、I型断裂、Ⅱ型断裂以及剪切摩擦实验，得到各自的荷载-位移曲线，计算弹性模量、峰值强度、断裂能、断裂位移等参数。材料参数如下表1～表3所示。

表1：实体单元力学参数

表2：裂隙单元力学参数

表3：接触对力学参数

S4：编制VUSDFLD子程序，自动删除煤集合单元，模拟工作面推进过程；

具体包括如下内容：

S4.1：在前处理程序中设置需要删除的单元集合

在步骤S2的数值模型中设置需要删除的单元集合，具体如下：

(1)首先确定数值模型中需要删除煤层单元的范围；

(2)对步骤S2中建立的单元集合“*Elset＝solid”、“*Elset＝coh”分别复制，然后通过减选方式框定需要删除的煤层单元；

(3)根据实体单元平均尺寸，以及工作面实际推进速度，确定每一次推进距离下需要删除单元的数量。然后通过减选方式，进一步确定每一次推进距离下需要删除的实体单元和裂隙单元，得到各自的单元编号；

S4.2：编制VUSDFLD子程序

第一，在编译器中，首先定义导入VUSDFLD子程序的参数，采用“subroutinevusdfld”子程序，导入主程序中的参数，包括单元编号、积分点坐标、状态变量参数、应力、位移、总时间、分析步时间、时间增量、方向等，并定义各参数为标量、矢量及张量；

第二，定义每一开挖步需要删除的单元集合；在定义好总时间totalTime、每推进一定距离所需要的时间deleteTime、单元号nblock、单元积分点坐标coordMp、方向direct的基础上，固定需要删除的起始单元积分点坐标coordMp(0)；

在此基础上，当deleteTime与totalTime的比值为c时，c为常数，定义所有满足该条件的单元为“集合1”；

当deleteTime与totalTime的比值为2c时，定义所有满足该条件的单元为“集合2”；依此类推；

第三，设置初始时刻所有单元的初始状态变量stateOld＝1，而后通过stateNew＝0命令，依次删除“集合1”、“集合2”等单元集合；由此实现单元删除、工作面推进的目标。

S4.3：将子程序导入主程序，模拟工作面推进过程

在Job模块下，新建立一个作业任务，在“通用”按钮下选择子程序文件，即可将步骤S4.2的VUSDFLD.for子程序导入主程序。

本实施例中，在数值模型中，根据实体单元平均尺寸，以及工作面实际推进速度，确定每一次工作面推进距离下需要删除单元的数量及其对应的边界位置。对步骤S2中设置的单元集合“*Elset＝solid”、“*Elset＝coh”分别复制。然后通过减选方式框定需要删除的煤层单元，并将需要删除的实体单元依次命名为solid-KW01，solid-KW02；将需要删除的裂隙单元依次命名为coh-KW01，coh-KW02。

在前处理的Materials模块中定义材料的状态变量参数Depvar＝1。在VUSDFLD子程序中，定义总时间totalTime＝10，每推进一定距离所需要的时间deleteTime＝5。第一次需要删除的实体单元编号nblock＝287～309，需要删除的裂隙单元编号nblock＝10103～10135；第二次需要删除的实体单元编号nblock＝462～490，需要删除的裂隙单元编号nblock＝10294～10326。当deleteTime与totalTime的比值为0.5时，stateNew(“1”,1)＝0，删除第一个“单元集合”；当deleteTime与totalTime的比值为1时，stateNew(“2”,1)＝0，删除第二个“单元集合”。由此实现采煤工作面向前推进的目的。

在Job模块下，新建立一个作业任务Job15，在“通用”按钮下选择子程序文件，将VUSDFLD.for子程序导入主程序。

S5：施加数值模型边界条件，输入材料力学参数，设置增量步大小，设置输出参量，包括，应力、应变、位移、单元节点力、裂隙单元的状态变量；

本实施例中，在Materials模块中输入材料参数(表1～表3)。在Assembly模块下对数值计算模型除上表面以外5个边界的法向位移进行约束，上表面为自由边界。而后根据采煤工作面的推进速度和模型尺寸，在Step模块下设定1个分析步，打开几何非线性，时间长度为10，时间增量设置为0.001；设置输出变量为应力、应变、摩擦位移、摩擦应力、状态变量STATUS等。通过自平衡方式确定模型的初始应力。

S6：通过连续-离散耦合模拟程序，求解模型，得到特厚煤层综放工作面底板裂隙的时空演化规律；

连续-离散耦合模拟程序的内容是：

在初始边界条件、和地应力条件下，通过有限元方法计算数值模型所有节点的节点力，通过步骤S4.2的VUSDFLD.for子程序删除单元后，采用有限元方法计算实体单元节点力，而后通过实体单元与裂隙单元的共享节点传递应力，利用有限元方法计算裂隙单元的位移量w；

采用离散元方法计算裂隙面的剪切摩擦效应，实体单元，即模拟完整煤岩块体、裂隙单元，即模拟非贯通裂隙、接触对，即模拟贯通裂隙的力学参数根据步骤S3确定；通过连续-离散耦合模拟程序，求解工作面推进过程中底板岩体应力、应变、位移等规律；

采动底板裂隙生成的数学判定准则为：根据上述连续-离散耦合模拟程序，计算裂隙单元的节点力、节点位移。通过下式计算任一时刻裂隙单元的断裂能：

其中，G

当有G

本实施例中，在Job模块下，对新建立的作业任务Job15，在Edit Job选项卡中设置并行计算核心数量为48核，提交计算。

S7：S7：采动底板导水通道演化分析

根据数值模拟结果，对超前支承压力、采空区底板卸压、特厚煤层综放动压这种循环加卸载条件下的底板导水通道发育规律进行分析，进而结合地下含水层分布情况，优化煤炭开采方案，实现绿色采煤。

本实施例中，数值计算得到特厚煤层综放开采工作面底板裂隙时空演化结果。通过连续-离散耦合程序，数值计算采动影响下底板各岩层裂隙单元的节点力、节点位移，并积分求解各自断裂能。以粉砂岩、中砂岩、煤、泥岩的I型和Ⅱ型断裂能实验结果为数学判定准则，对比数值求解断裂能与实验结果，判断是否删除裂隙单元。由此数值实现采动底板裂隙萌生、扩展、贯通的时空演化过程。

分析结果可知：

采动底板裂隙分布与演化规律，如图5黑色线条所示。从计算结果可以看出，模型中特厚煤层开采会引起其底板破坏。采空区范围内、靠近工作面位置处的底板，整体上按照塑性滑移线理论进行破裂；采空区范围内，距离工作面稍远位置处的底板，处于卸压状态；而在采空区范围内，远离工作面位置处，特厚煤层顶板垮落(图6和图7)。在这种循环加卸载(超前支承压力加载-采空区卸载-覆岩运动导致的动载)条件下，会进一步加剧底板损伤破裂。

在采空区之外且距离采煤面15m范围内，煤层出现张拉破坏，而底板出现拉/剪复合型破坏模式。且底板导水裂隙超前工作面煤壁25m产生，底板破坏深度达30m。相比较中厚煤层综采工作面条件，特厚煤层综放开采条件下底板破坏深度更大。

本实施例的特厚煤层综放工作面矿压显现的有限离散元模拟方法，与现有底板破坏数值方法相比，一方面实现了底板岩石在采动应力状态下从连续介质转化为离散介质的过程；另一方面，深刻反映了特厚煤层综放开采条件下，底板先后承受超前支承压力加载-采空区卸载-覆岩运动导致的动载这种静力-动力加卸载条件，以及由此引发的底板损伤破裂。为特厚煤层综放开采工作面底板导水裂隙演化提供了有力的数值模拟工具。

需要说明的是，以上的实施例是本发明较优的例子，本发明不限于上述实施例。任何熟悉本专业的技术人员，在不脱离本发明技术方案的范围内，作出添加、改进或修饰，以及等同变化，均属于本发明权利要求书所限定的保护范围内。