双系煤层采动导水裂隙演化规律的FDEM耦合模拟研究

李 浩，白海波，马立强，康志勤，**永，缪小成，武鹏飞，韦 婕

(1.太原理工大学 原位改性采矿教育部重点实验室，山西 太原 030024；
2.中国矿业大学 深部岩土力学与地下工程国家重点实验室，江苏 徐州 211116；
3. 中国矿业大学 矿业工程学院，江苏 徐州 211116；
4. 徐州中矿岩土科技股份有限公司，江苏 徐州 211116；

5. 晋能控股煤业集团，山西 大同 037003)

采动覆岩导水裂隙分布及演化规律的确定是开展水体下采煤、保水采煤、煤矿地下水库设计、断层防水煤岩柱留设、覆岩隔离注浆充填采煤设计等的重要基础[1-3]。长期以来，覆岩导水裂隙高度往往根据“经验公式”[4-5]来计算，一定程度上满足了煤矿安全开采的需求。然而，经验公式忽视了采动岩层连续→离散介质转化力学特性，导致在某些特定开采条件下的预计结果与实际偏差较大，甚至发生突水事故[6-7]。

建立适用于采动岩体连续→离散介质转化的耦合本构模型及相应的计算程序，是数值实现采动覆岩导水裂隙分布演化规律的重要手段。目前，学者广泛采用“连续介质假设”，通过有限元法(FEM)、有限差分法(FDM)，分别基于理想弹塑性模型[8]、应变软化塑性模型[9]、考虑围压对后继屈服函数影响的塑性模型[10]、弹塑性损伤模型[11]，以Mohr-Coulomb屈服准则[8-9]、应变能[12]、损伤变量[11,13]为判据，得到塑性区、损伤区高度为采高的6.6～27.5倍。此外，众多学者采用基于离散元方法(DEM)，将煤岩基质简化为矩形块、圆球形状并赋予弹性体或刚体力学属性，而基质间节理往往采用库伦滑移模型或法向/切向接触滑移模型[14-15]，以拉、剪应力为判据，实现采动覆岩破裂目的，根据牛顿第二定律计算得到覆岩裂隙带高度为采高的11.2～21.3倍。

但是，上述成果大多基于单一(系)煤层地质条件得到。更重要的是，在复杂采动应力状态下，完整覆岩往往由连续体转化为随机形状的离散体，在此过程中准脆性岩石出现拉-剪混合型韧性断裂、天然与新出现的结构面出现分离/压缩/剪切摩擦耦合响应。连续-离散转化过程深刻影响覆岩导水裂隙带分布与演化，仅基于“弹塑性、损伤、线弹性断裂力学、摩擦力学”的FEM或DEM等方法难以解决。

笔者紧密围绕采动覆岩连续-离散转化特性，依次建立岩石在复杂应力状态下的混合型韧性断裂本构、粗糙结构面压缩-剪切摩擦本构，编制相应的有限元-离散元(FDEM)计算程序。在验证本构模型合理性的基础上，模拟研究双系煤层采动导水裂隙演化规律。

受荷岩石的连续→离散转化过程包含3个显著的力学响应特征：① 拉、剪断裂为岩石的2种基本破坏类型，在复杂应力状态下往往形成混合型断裂模式；
② 岩石是准脆性材料，在断裂过程中表现出一定的韧性；
③ 完全断裂后，岩块间粗糙结构面可出现分离、压缩、剪切摩擦等复杂响应。以下建立岩石连续→离散介质转化的混合型韧性断裂-剪切摩擦(Mixed ductile fracture-shear friction，或MD-SF)本构模型，在此基础上编制有限元-离散元FDEM数值计算程序。

1.1 混合型韧性断裂本构

岩石韧性断裂力学响应通过有限元FEM方法计算。不同断裂模式下的断裂过程分为弹性变形和韧性断裂2个阶段。弹性段本构关系可表示为

(1)

一旦达到如下应力条件，岩石开始进入韧性断裂过程[16]：

(2)

为推导混合模式下材料完全断裂的判据，首先根据荷载-位移曲线(图1(c),(e))积分得到拉、剪及混合模式下的断裂能Gc,n,Gc,S,Gc,m，表达式为

(3)

式中，Gc,j为单位面积上黏聚力单元的断裂能，N/m;j为变量，表示n,S,m(即张拉型、剪切型、混合型)断裂模式。

其次，在复杂应力状态下，材料总是发生拉、剪混合型断裂模式，为确定材料完全断裂时拉、剪荷载的贡献，根据剪切型断裂能占总能量的比例，导出混合比ξ(无量纲)[16]表达式为

ξ=Gc,S/(Gc,n+Gc,S)

(4)

假设岩石在2个正交方向上的剪切性能相同，则材料在任意混合比ξ下完全断裂的判据[17]为

(5)

图2 拉、剪荷载下的荷载-位移曲线Fig.2 Load-displacement curves under tension and shear load

σc,n=(1-dc,n)Dc,nSc,n

(6)

(7)

(8)

(9)

由式(1)～(9)反映的岩石混合型韧性断裂力学性质可由图3表示,依次包括：① 弹性阶段应力-应变关系；
② 混合型韧性断裂起始时刻应力表达式；
③ 混合型断裂过程中应力-应变关系；
④ 混合模式下的完全断裂准则。

图3 混合型韧性断裂力学响应示意Fig.3 Schematic diagram of mechanical response of mixed ductile fracture

一旦断裂能满足式(5)时，岩石完全断裂。在数值计算中删除该断裂面上的黏聚力单元，从而消除其力学影响。

1.2 粗糙结构面压缩-剪切摩擦本构

岩石完全断裂后粗糙结构面的力学响应采用离散元方法(DEM)计算。在数值计算中，根据块体接触与相对运动趋势，将块体间关系划分为分离、压缩、剪切摩擦3类：① 当相邻实体单元(模拟完整岩块)的任意节点间距l>0，发生分离；
② 若l=0且2者间存在压应力时，块体相互挤压；
③ 若l=0且2者间存在沿结构面的剪应力，则出现剪切摩擦。其中，法向、剪切应力判别方法根据笔者前期研究进行[16]。

若相邻块体分离(l>0)，2者间不存在相互作用力，其运动符合牛顿第二定律[18]。

若相邻块体挤压，各要素本构关系为

(10)

式中，σn为压缩应力，MPa；
Nmax为结构面最大闭合量，mm，常参数，根据三维形貌扫描实验确定；
n为结构面闭合量，mm；
Dn为结构面法向模量，GPa，本文取值为其相邻岩块参数。

若l=0，以SARGIN[19]的研究成果为基础，进一步通过剪切实验曲线拟合得到粗糙结构面压剪摩擦本构关系为

(11)

式中，σS(σS,p)为(峰值)剪应力，MPa；
S(sp)为(峰值)剪切位移，mm；
DS为剪切刚度，GPa，3者均为常参数；
参数p由σS,p、Sp、残余剪应力σS,r和残余剪位移Sr确定，σS,p,σS,r满足如下强度准则[19]：

(12)

σS,r=σntanφr

(13)

其中，JRC为结构面粗糙度系数，无量纲；
σUCS为结构面抗压强度，MPa；
φb为结构面基本摩擦角，无量纲；
φr为结构面残余变形阶段的基本摩擦角，无量纲。需要强调的是，式(11)～(13)中的σp，DS，sp，p，σr，Sr等参数均为法向应力σn和结构面粗糙系数JRC的函数，其值根据图4实验获取。

图4 岩石断裂与剪切摩擦实验、数值模拟结果对比Fig.4 Comparison of experimental and numerical simulation results of rock fracture and shear friction

基于以上各式，建立混合型韧性断裂-结构面剪切摩擦(mixed ductile fracture-shear friction，MD-SF)本构方程组，编制相应的有限元-离散元(FDEM)数值程序，进而模拟岩石从连续介质→离散介质的转化过程(图1)。

岩石MD-SF本构方程的材料参数由断裂力学实验、结构面直剪实验获取。通过对比数值模拟结果与室内试验结果，验证本构方程合理性。

2.1 材料参数

岩石基质的力学参数根据文献获取[20]。张拉/剪切型断裂参数根据三点弯曲实验、加卸载条件下的贯穿剪切实验获取(表1)，基于各岩性试件的荷载-位移曲线，结合式(6)，(7)计算得到各自在拉、剪荷载下的损伤演化曲线(图1(d),(f))；
对于粗糙结构面剪切摩擦参数，通过如下方式获取：对岩心断口的三维形貌进行扫描，统计各岩性的JRC平均值，而后采用巴西劈裂方法制备粗糙结构面试件，选取与岩心断口相同JRC值的试件开展直剪实验。根据大同矿山4煤层埋深，确定法向应力σn=1，5，10 MPa。由此得到不同σn条件下粗糙结构面剪切摩擦参数(表1)。

表1 岩石力学参数Table 1 Rock mechanical parameters

2.2 本构方程验证

MD-SF本构方程的适用性通过偏置预制裂缝的三点弯曲、贯穿剪切、粗糙结构面直剪试验与模拟加以验证。

通过对比不同混合比ξ条件下的断裂力学实验与模拟结果，验证MD-SF本构方程中混合型韧性断裂理论。根据断裂力学实验(图1(a),(b))的试件尺寸、边界条件，建立图4(a)所示的数值计算模型。在模拟过程中，采用单调加载方式，位移加载速率为0.02 mm/min，加载时长10 min。模型的断面上布置零厚度的黏聚力单元，而在其他位置处布置实体单元。单元力学属性参见式(1)～(9)，力学参数由表1和图1(d),(f)确定。

对于直剪试验，以棱长100 mm立方体岩样为基础，采用巴西劈裂方法制备含粗糙结构面的试件，其JRC值为15(图4(b))。考虑到表1力学参数是根据法向应力σn=1,5,10 MPa得到的，为验证本构模型合理性，设置σn=3,8 MPa，固定下盘，上盘的剪切加载速率为0.02 mm/min，当残余应力稳定后停止加载。根据直剪试验条件建立相应的数值计算模型，结构面设置为平面，并在该处赋予法向、剪切摩擦接触属性(式(10)～(13))，其他位置处布置实体单元，力学参数由表1和图1确定。通过实验和FDEM计算程序(图5)，得到实验与模拟结果如图4所示。分析图4可得：

图5 FDEM数值计算流程Fig.5 Block diagram of finite element-discrete element numerical calculation

(2)将σn=1,5,10 MPa条件下实验所得力学参数(表1)代入MD-SF本构方程，通过图5所示计算流程得到σn=3,8 MPa下的剪应力-剪位移曲线(图4(d))。对比发现，实验与模拟结果吻合程度较好。具体来说，当σn从3 MPa增加至8 MPa时，实验(模拟)所得粗糙结构面峰值剪应力分别为4.82(4.77),10.68(10.94)MPa，峰值位移分别为0.89(0.93)，1.01(1.10)mm，残余剪应力分别为2.04(2.42)，5.63(6.02)MPa，残余剪位移分别为8.47(5)，7.13(5)mm。无论是模拟结果还是实验结果，均反映峰值强度、位移随σn的增加而增加，且模拟误差介于4%～8%。但模拟所得残余应力、残余位移与实验结果偏差较大，分别达到18%和69%。这是因为图4所示数值模拟目的在于判断MD-SF本构方程的合理性，考虑到计算成本，设置的计算时长只有250 min，残余位移最大为5 mm。尽管如此，MD-SF本构方程也适用于粗糙结构面的剪切摩擦过程。

(3)需要注意的是，图4(c)所示数值模拟所得曲线开头有非线性段。这是因为在数值模型位移加载前期设置5 min的线性加载阶段，加载速率逐渐从0增加至0.002 mm/min，避免形成冲击荷载。

3.1 工程概况与数值计算模型

大同矿区西北部侏罗系、石炭—二叠系可采煤层共计5层。侏罗系4个煤层开采深度约142～248 m，现已开采完毕。目前主采煤层为石炭—二叠系山4号煤，采深约420 m，采厚6 m。经过长时间、大规模地多煤层开采，侏罗系煤系地层已形成相互连通的导水裂隙，并在煤层采空区形成高达4×104m3采空区积水(图6)，严重威胁其下部山4号煤安全开采。

为研究大同矿区多煤层开采条件下采动覆岩运动及山4号煤的导水裂隙带发育规律，根据工程地质情况建立如图6所示的数值计算模型。考虑MD-SF本构方程的高度非线性，模型尺寸控制为x=870 m(长),y=450 m(高),z=1 m(厚)，在各自方向上分别约束其法向位移，且在x,z方向上设置侧压力系数分别为0.52，0.74，y方向设置自重载荷。工作面推进速度为8 m/d。覆岩中主要结构面用泰森多面体表示[21-22]，其间距根据如下方法确定：依据山4号煤顶板粉砂岩周期垮落步距确定其主要结构面的间距，再根据各岩性岩石与粉砂岩抗拉强度之比，确定粉砂岩、细砂岩、中砂岩、泥岩、高岭土砂岩、煤在x方向上的主要结构面平均间距(图6红色折线)分别为16.5，26.5，12.3，12.1，9.5，12.1，5.6 m，在y方向上的间距分别为10.3，20.0，10.6，5.0，4.7，6.0 m，若岩层厚度小于该值则取值岩层厚度，结构面初始宽度和JRC值见表1，而后赋予其“压剪摩擦(SF)”力学属性；
在此基础上，进一步在泰森多面体内部划分六面体实体单元，单元尺寸根据各岩性岩芯平均长度之比与计算成本综合确定，粉砂岩、细砂岩、中砂岩、泥岩、高岭土砂岩、煤的实体单元平均尺寸分别为4.9，4.1，3.8，5.0，3.6，5.6 m(图6黑色折线)，若岩层厚度小于该值则取值岩层厚度。实体单元赋予弹塑性本构模型[10]，而在单元边界嵌入零厚度黏聚力单元，并赋予其“韧性断裂-压剪摩擦(MD-SF)”力学属性，由此实现采动覆岩从连续体到离散体的转变。材料参数如表1和图1所示。

图6 采空区积水分布及采动覆岩运动数值计算模型Fig.6 Water distribution in goaf and the numerical calculation model of mining overburden movement

3.2 模拟结果

通过上述简化模型，得到采动覆岩导水裂隙带数值计算结果(图7)。由上述数值模拟结果可得：

图7 不同工作面推进距离下岩石连续-离散转化、覆岩裂隙分布及Mises应力场特征Fig.7 Characteristics of rock continuous-discrete transformation,overburden fracture distribution and Mises stress field under different mining face advancing distance

(1)在工作面推进方向上，覆岩导水裂隙介于工作面煤壁前方44.7 m与煤壁后方12.4 m之间。岩体脆-韧性力学特征、超前支承压力、顶板天然结构面位置的叠加效应是产生上述现象的主要原因。山4号煤开采过程中，超前支承压力范围约为49 m，这会导致该范围内脆性砂岩的结构面承受较大压剪应力；
随着工作面向前推进，顶板岩体垮落，其附近岩体失去侧向支撑而旋转，导致在脆性明显的砂岩内较大宽度的超前导水裂隙形成；
反之，若工作面煤壁刚好经过前一条导水裂隙，且超前支承压力不足以使后一条结构面出现法向位移时，则出现滞后型导水裂隙。与此同时，2条天然结构面之间的完整岩体内部出现复杂的拉、剪应力，导致完整岩体主要产生混合型断裂裂缝，并最终由连续体转化为离散体，使得覆岩导水裂隙进一步复杂化。

(2)采动覆岩导水裂隙高度。在山4号煤覆岩内且与煤层间距45，82，125 m处布置3条导水裂隙宽度测线(图6)，统计裂隙宽度在1mm以上的导水裂隙总宽度随工作面推进距离演化规律如图8所示。

图8 覆岩导水裂隙总宽度随工作面推进距离曲线Fig.8 Curves of total width of water conducting fissure in overburden with mining advancing distance

图8中上、中测线位于高岭土砂岩层中，在工作面推进过程中，2个测线处导水裂隙总宽度呈现先增加后减小的趋势，即从0增加至0.476 m(或1.053 m)再降低至0.143 m(或0.520 m)。以图7(b),(c)为例解释该现象的原因：当工作面推进至123 m时出现周期来压，顶板砂岩层完全断裂并发生约4°的转动，裂隙宽度迅速增加至0.173 m；
而随着工作面推进至184 m，采动应力逐渐稳定，该岩层回转约2°，使得裂隙宽度降低至0.134 m。以上数据表明，煤层开采过程中更容易出现顶板突水事故。

由图7,8可知，高岭土砂岩中的采动覆岩裂隙宽度、密度远远小于其下部砂岩层。这不仅是因为砂岩层受采动应力较大，更重要的是，高岭土砂岩具有更强的韧性破坏特征，在断裂过程中大量能量转化为塑性功，而高岭土砂岩的裂缝表面能(即弹性能)是砂岩的63.57%，这意味着采动裂隙在脆性较为明显的砂岩内更容易扩展。此外，砂岩较为显著的脆性断裂响应，导致采动裂隙快速发育至砂岩层顶界，岩块偏转并与周围未破断岩体形成最宽达3.347 m的宽大裂隙，而随着采动应力降低岩块回转，裂隙宽度又会逐渐降低，并稳定在1.726 m左右，形成波动升降曲线(图8蓝色线)。

综上，山4号煤的采动覆岩导水裂隙带高度不超过125 m，远大于经验公式所得74.1 m的结果。

为验证覆岩导水裂隙带模拟结果，进一步开展井下注水试验。钻孔距离开切眼280 m，仰角65°，钻孔斜长135 m(垂高120 m)，待工作面开采完毕后采用分段注水方法得到水流漏失量。由图9可知，当钻孔斜长从7.6 m增加至23.4 m时，水流漏失量从2.2 L/min快速增加至26.5 L/min，意味着钻孔进入导水裂隙带；
直至钻孔斜长达到117.6～123.9 m，水流漏失量为4.2～2.7 L/min，达到稳定。实测结果表明，山4号煤采动覆岩导水裂隙带垂高约为112.3 m，从而验证上述模拟结果的合理性。

图9 注水试验结果Fig.9 Water injection test results

(1)建立混合型韧性断裂-剪切摩擦本构模型及相应的FDEM数值算法，采用实体单元、黏聚力单元和接触对表示岩石基质、非贯通裂隙及贯通结构面，计算准脆性岩石弹塑性变形、混合型韧性断裂及粗糙结构面分离/压缩/剪切摩擦力学响应，数值实现岩体从连续介质到离散介质的转化过程。

(2)相比较实验结果，数值模拟所得不同拉、剪断裂能混合比下岩石断裂峰值荷载、峰值位移、断裂能误差为0.5%～3.0%，2.0%～11.0%，4.0%～16.0%；
模拟所得不同法向应力下粗糙结构面剪切摩擦峰值荷载误差为4%～8%，由此验证混合型韧性断裂-剪切摩擦本构模型的合理性。

(3)模拟条件下，开采过程中覆岩导水裂隙总宽度是覆岩运动稳定后的2.26～7.11倍，表明煤层开采过程中更容易出现顶板突水事故。

(4)高岭土砂岩显著的韧性破坏特征有效控制其内导水裂隙发育高度，但在双系煤层重复开采扰动下，覆岩导水裂隙发育高度约125 m，裂采比为20.8，远超经验公式所得结果。

推荐访问:裂隙 煤层 耦合