采动应力下煤体渗透率模型构建及研究进展

周宏伟，荣腾龙，牟瑞勇，王路军，任伟光

(1.中国矿业大学(北京) 能源与矿业学院，北京 100083; 2.中国矿业大学(北京) 力学与建筑工程学院，北京 100083; 3.中国矿业大学(北京) 煤炭资源与安全开采国家重点实验室，北京 100083)

煤炭是我国的基础能源，也是重要的工业原料，到2030年，煤炭科学开采预测产能总量达30～35亿t[1]。所以，未来一段时间内煤炭生产仍将持续进行，煤矿的安全生产仍需加强保障。煤炭形成过程中伴生有大量的瓦斯气体，瓦斯既是瓦斯爆炸、煤与瓦斯突出灾害的诱因，又是一种高效的清洁能源[2-4]。掌握煤层中瓦斯流动规律是减少瓦斯事故、扩大瓦斯采集利用的前提，而煤体的渗透率是决定瓦斯流动规律的关键参数[5]，所以对煤体渗透率进行研究是十分必要的。我国大部分矿区的煤层渗透率较低，需通过采动卸压、水力压裂、深孔爆破等手段进行人工增透[6-8]，此外煤与瓦斯突出事故多发于掘进工作面[9-11]。因此，建立工作面前方采动应力下的煤体渗透率模型对提高瓦斯抽采和减少瓦斯事故具有重要意义。

以往的煤体渗透率模型较多以煤层气开采为背景[12-14]，同时基于单轴应变条件和上覆载荷不变的假设建立，而在煤矿采掘活动中，煤层的受力特点有别于煤层气开采。工作面推进过程中采煤工作面前方煤体经受采动应力路径的变化过程，其竖直方向的支承压力逐渐升高后下降、水平方向的应力持续降低，煤层受力呈现典型的三向应力变化特征(图1(a))。随着支承压力的升高以及人工增透活动的进行，煤体内部继发出现损伤和破裂现象，煤层内裂隙扩展贯通造成瓦斯压力下降，同时伴随瓦斯解吸扩散的发生，最终煤体渗透率急剧增长。采动应力决定了煤体在由原岩应力状态渐近为侧面临空过程中的受力状态、物理力学特征、温度梯度及渗透率变化有别于煤层气的抽采过程，以往的煤体渗透率模型并不能完全适用于采动应力条件。

笔者首先概述了煤体渗透率经典模型，然后结合采动过程中煤体内的力学变化机制及渗透率的控制因素提出了采动应力下煤体渗透模型构建中的3个关键问题，即采动煤体各向异性特征、采动煤体损伤破裂特征和煤体吸附解吸特征，并将这3个方面的国内外研究成果进行了总结和分析。在此基础上，将各向异性特征的煤体渗透率模型归纳为由有效应力变化表征和几何参数变化表征的两类，又将损伤破裂特征的煤体渗透率模型划分为本构方程中含损伤变量和渗透率表达式中含损伤变量的两类。最后，对采动应力下煤体渗透率模型的研究进展进行了介绍并对当前采动应力下煤体渗透率模型研究的不足和未来的发展方向进行了探讨。

1 煤体渗透率经典模型概述

为了对煤体的渗透率演化过程进行定量描述，国内外学者建立了大量的渗透率模型，这些研究成果主要集中在煤层气开采方面，其中一些模型具有较高的代表性，这些经典的模型成果在煤体渗透率模型的研究领域中发挥了重要作用。

GRAY[15]认为煤储层的渗透率与垂直于割理结构的有效水平应力大小有直接关系，较早建立了煤体渗透率演化模型。SEIDLE等[16]将储层煤体简化为“捆绑的火柴棍”几何结构，得到了煤体渗透率与有效静水压力的关系，从理论角度推导出渗透率变化的指数函数形式。PALMER和MANSOORI[17]假设储层为单轴应变条件，考虑煤层割理压缩性的变化，基于渗透率k与孔隙率φ的立方关系，结合非等温多孔弹性理论建立了煤层渗透率模型。SHI和DURUCAN[18]从各向同性的弹性本构方程出发，将煤储层视为连续介质，基于单轴应变的假设条件求出煤储层水平方向有效应力的变化，建立了指数函数形式的渗透率模型。CUI和BUSTIN[19]考虑整体平均有效应力对煤体孔隙度和渗透率的影响，利用多孔弹性本构方程进行推导得到解吸过程中煤体渗透率模型。ZHANG等[20]将煤体视为固体颗粒和孔隙空间组成，假设煤的孔隙空间同样存在吸附应变，建立了一个广义孔隙度的渗透率模型。LIU等[21]认为煤体渗透率由基质渗透率和裂隙渗透率2部分组成，考虑基质和裂隙的相互作用得到合成的渗透率表达式。经典的煤体渗透率模型汇总见表1。

表1煤体渗透率经典模型
Table1Classicalpermeabilitymodelsofcoals

文献来源模型表达式参数释义GRAY[15]σeh-σeh0=-ν1-ν(p-p0)+E(1-ν)ΔεsΔpsΔpsk=1.013×10-0.31σehσeh为水平有效应力;ν为泊松比;E为弹性模量;p为气体压力;Δεs/Δps为单位等效吸附压力改变下的煤体应变;Δps为等效吸附压力变化增量;k为煤的渗透率;下标0为初始状态SEIDLE等[16]kk0=exp[-3Cf(σehy-σehy0)]Cf为割理压缩系数;σehy为有效的静水压力PALMER和MANSOORI[17]φφ0=1+cmφ0(p-p0)+clφ0KM-1 Bp1+Bp-Bp01+Bp0 kk0=φφ0 3M=E(1-ν)(1+ν)(1-2ν),K=E3(1-2ν),cm=1M-KM+G-1 γs,其中γs为颗粒压缩系数;φ为孔隙率;cl为Langmuir吸附应变常数;B为Langmuir吸附平衡常数;G为比例系数SHI和DURUCAN[18]σeh-σeh0=-ν1-ν(p-p0)+E3(1-ν)εLpp+pL-p0p0+pL kk0=exp[-3Cf(σeh-σeh0)]εL为Langmuir吸附体积应变常数;pL为Langmuir吸附压力常数CUI和BUSTIN[19]kk0=exp-3Kp[(σm-σm0)-(p-p0)] Kp为孔隙体积模量;σm为平均法向应力ZHANG等[20]kk0=11+S(1+S0)+αbφ0(S-S0) 〛 3S=εV+(p/Ks)-εs,S0=(p0/Ks)-εLp0/(p0+pL),αb=1-K/Ks,Ks为煤固体颗粒体积模量;K为煤的体积模量;εV为煤的总体积应变;εs为煤的吸附应变LIU等[21]kk0=km0km0+kf01+Rmpmφm0K 3+kf0km0+kf01+1-Rmφf0(ΔεV-Δεs) 〛3km为基质渗透率;kf为裂隙渗透率;φm为基质孔隙率;φf为裂隙孔隙率;Rm为整体弹性模量损失率;pm为基质上的气体压力

2 采动应力下煤体渗透率模型构建的关键问题

2.1 采动煤体各向异性特征

煤属于沉积岩，在沉积方向具有明显的层理结构，而在垂直层理的方向又存在空间狭窄的裂隙结构，这种结构被称为割理结构。割理可分为面割理和端割理结构，通常将两种割理结构简化为相互正交的形式。因此，煤体具有显著的各向异性特征(图1(b))，进而煤体在各个方向上存在不同的渗透率变化趋势。除此之外，采动应力下煤体渗透率模型的构建必须把握三向受力的特点，建立可表征竖向加卸载、水平卸载过程的渗透率理论模型。因此，无论是煤体内部结构特征还是外部受力条件均要求使用各向异性的渗透率模型对采动应力下煤体的渗透率演化过程进行描述。笔者根据渗透率模型的建模角度将现有的各向异性煤体渗透率模型成果归纳为有效应力变化和几何参数变化两类进行介绍。各向异性特征的煤体渗透率模型汇总见表2。

表2各向异性渗透率模型
Table2Anisotropicpermeabilitymodels

文献来源模型表达式参数释义CHEN等[22]ki=ki0exp[3Cfi(σei-σei0)] i=x,yCfi为不同方向的割理压缩系数;σei为不同方向的有效应力LIU等[23]ki=ki0exp[-3Cfi(σem-σem0)] i=x,y,zCfi为不同方向的平均割理压缩系数;σem为有效的平均法向应力ZHANG等[24]ki=ki0exp[-αiCpσ(σi-σi0)] i=x,y,zCpσ为孔隙压缩系数;αi为幂函数系数;σi为不同方向的应力PAN和CONNELL[25]ki=ki0exp[3Cfi(σei-σei0)] i=x,y—GU和CHALATURNYK[26]kiki0=1+ajbjΔεfj 31-Δεfj i=x,yaj为j方向的基质单元宽度;bj为j方向的裂隙单元宽度;Δεfj为j方向的裂隙应变增量马天然等[27]ki=∑3j=1(bj0+Δbj)312(aj0+Δaj) i=x,y,z—LIU等[28]kiki0=∑i≠j121+3(1-Rm)φf0(Δεtj-13Δεs) 〛Δεtj为j方向煤体的总应变YANG等[29-31]kxkx0=121+3(1-R2)φ0 〛(Δεt2-13Δεs) 3+ 121+3(1-R3)φ0 〛(Δεt3-13Δεs) 3kyky0=121+3(1-R1)φ0 〛(Δεt1-13Δεs) 3+ 121+3(1-R3)φ0 〛(Δεt3-13Δεs) 3kzkz0=121+3(1-R1)φ0 〛(Δεt1-13Δεs) 3+ 121+3(1-R2)φ0 〛(Δεt2-13Δεs) 3Ri为i方向的模量损失率;Δεti为i方向的煤体总应变WANG等[32]kiki0=1+2(1-Rm)φf0Δεej 〛3 i,j=x,y,z;i≠jΔεej为j方向的煤体有效体积应变WU等[33]kmkm0=1(1+S)φm0[(1+S0)φm0+αb(S-S0)] 3kiki0=∑i≠j121-1φf0+3KfK13αTΔT+13Δεs-1KΔσei 〛3Kf为裂隙刚度;αT为热膨胀系数;T为温度WANG等[34]ki=ki0 1+1φi0Δσej-vbjkΔσek-vbijΔσeiEbj-Δσej-vmjkΔσek-vmijΔσeiEmj+Δσek-vbkiΔσei-vbjkΔσejEbk-Δσek-vmkiΔσei-vmjkΔσeiEmk-εLjFljppLj+p-Fljp0pLj+p0 -εLkFlkppLk+p-Flk0p0pLk+p0 〛 3i≠j≠kEbi为i方向的煤体弹性模量;Emi为i方向的基质弹性模量;εLi为i方向的Langmuir吸附体积应变常数;pLi为i方向的Langmuir吸附压力常数;Fli为i方向的内膨胀应变系数

2.1.1 有效应力变化表征的煤体渗透率模型

有效应力变化表征的煤体渗透率模型都以煤体所受有效应力为自变量，在建模过程中较多将煤基质和裂隙视为一个整体，利用弹性本构方程对煤体在坐标轴3个方向的应力和应变关系进行推导，建立渗透率和有效应力之间的关系。此类模型[22-25]均可以归纳为以下形式

(1)

(2)

2.1.2 几何参数变化表征的煤体渗透率模型

几何参数变化表征的煤体渗透率模型大多以煤体内孔裂隙尺寸、裂隙应变和孔隙度为自变量，建模过程基本是将煤体划分为基质和裂隙两部分，通过计算应力、吸附作用下裂隙尺寸或孔隙度变化，建立煤体裂隙开度、孔隙度和渗透率的关系。为了考虑煤体各向异性的结构特征(图1(b))，这些渗透率模型[26-34]构建过程中大多将煤体抽象为立方体几何结构模型(图1(c))。

一些学者将气体在煤层裂隙中的流动简化为平板间流体流动过程，通过POISEUILLE定律得到渗透率与裂隙开度之间的立方关系(式(3))，并以此开始推导,即

(3)

式中，a，b分别为基质单元和裂隙单元的宽度。

GU和CHALATURNYK[26]将煤体简化为火柴棍模型，而马天然等[27]将煤体简化为立方体模型，分别对式(4)进行变换直接得到了渗透率和裂隙开度的3次方关系式。LIU等[28]将煤体抽象为立方体几何结构，认为煤体为各向同性的线弹性材料，提出了模量损失率来表征各方向弹性模量的变化过程，建立了不同边界条件下的各向异性渗透率模型。YANG等[29-31]认为工程分析中应忽略基质渗透率的影响，借此对LIU等[28]模型加以修正，构建出横观各向同性的渗透率模型。WANG等[32]从式(4)出发对压缩过程中裂隙孔径的变化进行推导，得出一个与LIU等[28]模型形式类似的表达式，同时定义了基质与整体单元应变比的新参数，建立了一个2维各向异性的渗透率模型。

另外，一些学者定量表示出三向有效应力变化过程中煤体内部孔隙度的变化，并将孔隙度变化结果代入渗透率和孔隙度的立方关系[17](式(4))建立了相应的各向异性渗透率模型。WU等[33]将煤体的物理力学性质视为各向同性，假设坐标轴3个方向下热膨胀和吸附作用产生的煤体应变相等而有效应力不相等，建立了双孔隙度的煤体各向异性渗透率模型。WANG等[34]基于立方体模型提出了内膨胀应变系数以对吸附过程中煤体基质的内膨胀量进行修正，建立了一个以三向有效应力和气体压力为自变量的渗透率模型。

(4)

式中，k为渗透率;φ为孔隙率。

从建模过程上来说，有效应力变化和几何参数变化的两个建模角度可认为是互逆的相反过程。有效应力的角度是将煤体所受有效应力通过孔隙度、应变、模量、泊松比、吸附常数等参数综合表示出来，而几何参数的角度是将裂隙开度或孔隙度用应力、模量、泊松比、吸附常数等参数进行表达，进而各自代入相应的框架基础公式。此外，根据各自建模所依据的框架基础，有效应力角度建模结果基本为指数型函数、几何参数角度建模的结果多数为3次方的幂函数。

渗透率模型中如果能将煤体的边界条件和结构参数均考虑为各向异性且持续变化的特征是最为理想和贴近实际的情况，然而受目前试验条件和测试手段所限，大量学者也只能将模型中的部分参数做出近似和假设处理。从上述各向异性的模型成果中可以看出，外部边界的三向应力和三向应变参数较易跟踪和获取，所以三向应力的特征在模型中容易实现[22-25]。结构异性参数大多较难准确把握，比如3个方向的割理压缩系数、弹性模量和泊松比等，因此一些模型中将结构参数视为各向同性[26-33]。

2.2 采动煤体损伤破裂特征

受开采扰动的影响，支承压力的空间分布具有分区化特征，煤体由原岩应力逐渐升高至应力峰值后迅速卸压(图1(a))，在此过程中煤体由原生初始状态渐进为损伤破裂状态，伴随煤体内微裂纹的扩展贯通，渗透率将会迅速增大，可达初始值的数十倍，乃至百倍[35-36]。因此，对于采场前方煤体必须考虑开采扰动引起的损伤破裂对渗透率的影响。

目前，国内外学者已经开始探究损伤破裂影响下渗透率的变化规律。开采扰动下煤体中瓦斯流动及渗透率的变化是一个典型的气固耦合问题。一种解决此问题的手段就是将应力场和渗流场的控制方程进行耦合求解。根据损伤力学理论可以定义多种形式的损伤变量，损伤变量对渗透率的控制作用通过耦合方程的形式进行传递。这种调控方式下，损伤变量处于固体力学的本构方程之中，而渗透率表达式内并不含损伤变量。另一种处理此问题的方式就是将定义出的损伤变量直接嵌入渗透率模型的表达式之中。笔者依据损伤变量所处的位置不同将损伤特征的实现方式归纳为两种:一种为损伤变量存在于本构方程中，另一种为损伤变量存在于渗透率表达式中。考虑损伤破裂特征的渗透率模型成果汇总见表3。

表3考虑损伤破裂特征的渗透率模型
Table3Damage-basedpermeabilitymodels

文献来源模型表达式参数释义YANG等[48]k=k0e-β1(σ1-α1p) D=0ξk0e-β1(σ1-α1p) D>0 k=k0e-β1(σ3-α1p) D=0ξk0e-β1(σ3-α1p) 0

2.2.1 本构方程中含损伤变量的煤体渗透率模型

一些学者基于宏观唯象的连续损伤力学理论建立了煤体损伤本构并将其代入固体力学控制方程，与瓦斯渗流方程联立进行求解。基于LEMAITRE[37]提出的应变等效假设(式(5))，魏明尧等[38]建立了弹塑性损伤本构;HU等[39]定义了双重孔隙结构下煤体的等效弹性损伤柔度张量，构建了双重孔隙介质损伤的本构模型;WANG等[40-41]将HU[39]定义的等效弹性损伤柔度张量进行简化后应用在单一孔隙率的渗流方程之中,即

(5)

另外，一部分学者通过细观损伤力学的方法建立了煤体的损伤本构，将煤体视为大量代表性体积单元(RVE)组成，通过研究每个RVE的力学破坏特征来反映煤体内损伤演化过程。ZHU和TANG[42-43]考虑岩石的非均质性，假设岩石参数符合Weibull分布(式(6))，将最大拉应力和Mohr-Coulomb准则作为判断煤体RVE破坏的准则(式(7))，建立了煤体的弹性损伤本构。沿用这种研究思路，一些学者建立了不同的气固耦合模型用以探究煤层中损伤-应力-渗流的问题，间接搭建了损伤变量和渗透率之间的关系[44-52]。除了根据渗透率与孔隙度立方关系(式(4))推导渗透率表达式的常见方法[44-45]之外，YANG等[48]利用分段函数依次建立了煤体单元在压缩和拉伸条件下的渗透率函数式;刘力源等[51-52]将煤体视为双重孔隙介质，结合煤体吸附CO2出现的损伤演化过程分别得到煤基质和裂隙的渗透率计算式。

(6)

式中，C1，C2为两个损伤破裂阈值函数;u为单元参数;u0为单元参数的平均值;m为函数的形状参数;σ1,σ3分别为单元的最大主应力和最小主应力;ft0为单轴抗拉强度;fc0为单轴抗压强度;φ′为内摩擦角。

2.2.2 渗透率表达式中含损伤变量的煤体渗透率模型

通过理论推导构建损伤变量并将其嵌入渗透率表达式的研究主要限于连续损伤力学的框架内进行，大多数学者基于应变等效假设展开研究，而损伤变量的构建方法多样。XUE等[53]通过孔隙度和有效应力定义了损伤变量;ZHANG等[54]采用煤体的体积应变、吸附应变和初始孔隙率构建了损伤变量;LI等[55]以轴向应变构建了损伤变量。另外，也有学者采用回归分析的方法，通过拟合系数来确定渗透率和不同损伤程度下应力状态之间的关系。CHEN等[56]为了描述煤体从弹性状态到峰后破坏的渗透率变化全过程，提出了修正系数L，通过回归函数得到有效偏应力与修正系数L的“S”形关系。ZHANG等[57]基于达西定律和裂隙的本构方程，通过拟牛顿算法的回归分析确定了以应力表示的渗透率函数。为了搭建渗透率和损伤变量之间的“桥梁”，另外一些学者立足于试验结果反映出的变化规律，在前人建立的渗透率模型基础之上进行改造，直接添加反映损伤影响的作用项或修正系数对渗透率变化过程进行推演，也取得了一些研究成果。ZHU等[58]在渗透率和孔隙度的立方关系(式(4))上进行改造，添加了指数形式的损伤作用项;ZHENG等[59]以CUI和BUSTIN[19]模型为基础，在该模型中添加了损伤作用项，假设了渗透率和损伤变量的指数关系;XUE等[60]通过室内试验的渗透率和声发射数据反分析提出在峰后阶段引入一个作用项来反映损伤影响，并假设峰后阶段渗透率和损伤变量呈线性关系;XIE等[61]假设损伤作用下渗透率的变化符合式(8)的关系,即

(8)

式中，ζ为与材料有关的系数。

可以看出，考虑损伤破裂特征的渗透率模型研究在近十年才逐渐展开，该领域的研究尚处于起始阶段，模型中损伤效应的研究深度尚浅，损伤变量的构建方式也较为分散。在模型适用性上，损伤变量存在本构方程中的实现方式可以适用较多边界条件情况，具有更广的适用范围，但是该方式下求解煤体渗透率需要依次计算应力场和渗流场控制方程，计算工作量较大，所以该方式更适合采用数值计算。损伤变量存在渗透率表达式的实现方式使得研究渗透率与各因素之间的数量关系更加便捷，但其大多基于已有模型成果进行改进，受原模型构建过程中的假设条件限制，诸如单轴应变条件、上覆载荷不变等。损伤变量所处不同位置的两种处理方式都可以达到在渗透率计算过程中考虑损伤破裂影响的目的。究其本质而言，损伤变量存在于本构方程的情况可视为损伤变量存在渗透率表达式情况的隐式函数形式。

2.3 煤体吸附解吸特征的表征方法

不同于其他气体在岩石中的渗透行为，瓦斯在煤层中运移时存在吸附解吸的特殊现象。煤层采掘过程中，煤体损伤破裂伴随着瓦斯渗流通道的形成。与此同时，开采扰动打破了煤层内原有的瓦斯压力平衡状态造成瓦斯压力持续下降，煤体内吸附态瓦斯释放引起解吸现象的发生。吸附解吸作用直接影响煤基质的胀缩变形，在渗透率模型构建的过程中，大多通过吸附应变这一中间变量间接反映煤体所受有效应力的大小或煤体内部渗流通道几何尺寸的变化[62]，因此关于吸附应变的计算方法就成为了研究吸附解吸影响下渗透率变化的焦点。

在早期的煤体渗透率模型中，通常将吸附应变与气体压力的增量视为线性关系进行处理。GRAY[15]通过提出等效吸附压力的概念建立了吸附应变与气体压力之间的线性关系。LEVINE[63]利用基质收缩系数确定了基质吸附应变和气体压力增量之间的线性函数。国内外学者普遍将煤对瓦斯的吸附认为是物理吸附过程且为单分子层吸附。基于此，Langmuir等温吸附模型[64]可以较好的用于煤对瓦斯气体吸附量的计算(式(9))。HARPALANI和CHEN[65-66]认为在气体吸附过程中吸附应变和气体吸附量成正比，采用Langmuir理论公式计算瓦斯气体吸附量，建立了煤体吸附应变和瓦斯气体吸附量的线性关系(式(10))。鉴于Langmuir等温吸附模型可以较好的描述气体压力改变下煤对瓦斯的吸附气体量，LEVINE[63]对Langmuir模型作出进一步改进，采用一个和Langmuir吸附模型具有相同数学函数形式的方程(式(11))来描述煤体吸附应变。

εs=αsVs(10)

式中，Vs为吸附的瓦斯气体含量;VL为Langmuir体积常数;αs为煤的吸附膨胀系数;εL为Langmuir吸附体积应变常数;pL为Langmuir吸附压力常数。

煤体对瓦斯吸附解吸的内在机理可通过热力学理论进行探究，所以除了采用上述Langmuir形式的方程对吸附应变进行表征之外，一些学者通过吸附热力学方法建立了吸附应变的计算式[67-72]，详见表4。

表4吸附热力学角度建立的吸附应变表达式
Table4Expressionsofadsorptionstrainbasedonadsorptionthermodynamics

文献来源吸附应变表达式参数释义PAN和CONNELL[67]εs=RTAln(1+Bp)ρmEmf(x,νm)-pEm(1-2νm)R为气体常数;T为温度;A为Langmuir吸附常数;ρm为煤基质密度;Em为基质的弹性模量;νm为基质的泊松比LIU和HARPALANI[68]εl=3AρmRTEAVg∫p0B1+Bpdpεl为吸附/解吸引起的线性应变;EA为吸附解吸的膨胀模量;Vg为气体摩尔体积LIU等[69]εs=nsatRTVcEc[ln(1+Bpt)-ln(1+Bpr)]nsat为饱和吸附摩尔数;Vc为煤的体积;Ec为煤的弹性模量;pr为起始气体压力;pt为终止气体压力VERMOREL和PIJAUDIER-CABOT[70]εs-εs0=∫pp0dpK1-KKs (1-χ)-1-1 〛χ为约束度HOL等[71]εsij=Vm0ρ(Aij+δij)3Cσεsij为各个方向的吸附应变;Vm0为吸附相的偏摩尔体积;ρ为煤的密度;Aij为各向异性膨胀张量;δij为克罗内克符号;Cσ为吸附分子的摩尔数

采用以气体压力增量为自变量的线性函数对煤吸附应变进行描述是最早的研究方法，其优点是瓦斯气体压力的测定较为简单，但是线性函数却有夸大瓦斯压力对吸附应变影响的缺点。通过对比试验结果发现Langmuir等温吸附理论可以很好的描述煤对瓦斯的吸附解吸过程，上述Langmuir形式的吸附应变表达式具有一定的理论基础且方程形式较为简洁，表达式中的VL,εL,pL等参数也较易测试，所以式(10)[73-76]和式(11)[77-80]被后来学者大量效仿使用。通过热力学方法建立的吸附应变大多以吸附热力学为基础进行推导，由此得出的吸附应变表达式具有更为严密的理论基础，然而此类方法的模型参数较多且部分参数难以准确测定，所以此类吸附应变表达式的实用性有所限制。

3 采动应力下煤体渗透率模型研究进展

现有的煤体渗透率模型中同时考虑各向异性、损伤破裂、吸附解吸等特征的渗透率模型还较少。近些年，部分学者已经开始关注这一问题，并建立了相应的采动应力下煤体渗透率模型。笔者将这些模型成果归纳为有效应力型、几何参数型和系数拟合型的3种进行介绍。

3.1 有效应力型的煤体渗透率模型

有效应力型的煤体渗透率模型主要以有效应力、气体压力、弹性模量和泊松比等力学参数为基本变量对渗透率的变化进行表征。其中，ZHANG等[81]基于立方体模型结构(图1(c))，运用损伤和断裂力学理论提出了一个可考虑采动应力引起微裂纹扩展、吸附解吸造成煤体变形的有效应力型渗透率模型(式(12))，并根据谢和平等[82]提出的3种典型开采方式下采动应力路径(图2)，对采动煤体渗透率变化规律进行了拟合分析。

图2(a)～(c)为不同开采方式下煤体渗透率的试验和拟合结果。可以看出在3种开采方式下，随着采动卸压过程的进行，煤体内微裂隙不断萌生和扩展，损伤积累不断增大，渗透率的实测和拟合结果持续升高。无煤柱和放顶煤开采方式的渗透率在加卸载途中出现陡然上升的情况，但拟合结果却没有出现陡然升高的特征。不过整体来看，式(12)基本可以反映出煤体渗透率随着采动卸压过程的进行而持续升高的变化规律。

3.2 几何参数型的煤体渗透率模型

几何参数型的煤体渗透率模型主要以孔隙率和煤体应变为基本参数来构建渗透率模型。ZHENG等[83]将LIU等[28]提出的煤体各向异性渗透率模型和ZHU等[58]定义的损伤破裂特征修正项进行联结。通过细观损伤力学的方法构建损伤变量以对渗透率进行修正，合成出可表征采动应力下煤体损伤的几何参数型渗透率表达式:

图2 不同开采方式下应力路径及煤体渗透率演化[81-82]Fig.2 Stress paths and permeability of coal under different mining layouts[81-82]

式中，εj为j方向的应变;εt0,εc0分别为最大拉伸主应变和最大压缩主应变;ε1,ε3分别为最大主应变和最小主应变。

其后，利用式(13)分别对考虑和不考虑损伤效应的保护层底板下煤岩体渗透率进行了模拟分析。由图3可以看出距离保护层底板下方80 m处的岩体渗透率为初始状态;在与底板距离缩进过程中煤岩体渗透率逐渐升高;在距离底板12 m之外，考虑损伤与否的渗透率结果相同;而在12 m范围之内，考虑损伤效应的煤岩体渗透率比率出现急剧上升，从初始渗透率比率的180倍陡增至680倍左右，而不考虑损伤的煤岩体渗透率比率仅升高至200倍左右。

图3 保护层底板下煤体渗透率分布[83]Fig.3 Permeability distribution of the coal under the first-mined seam floor[83]

3.3 系数拟合型的煤体渗透率模型

系数拟合型的煤体渗透率模型主要通过试验数据拟合的手段来确定一系列待定系数的变化与不同区段渗透率的关系。XUE等[84]根据煤体三轴压缩试验过程中应力-应变-渗透率关系的变化特征(图4)，将全应力-应变过程划分为3阶段:弹性阶段、塑性软化阶段、残余状态阶段。同时认为弹性阶段没有微裂纹产生，塑性软化阶段萌生大量裂纹，进而结合巷道开挖后围岩损伤的3区特征，将围岩的渐进损伤区与全应力-应变过程3阶段进行了一一对应:弹性阶段对应无损区、塑形软化阶段对应开挖损伤区、残余状态阶段对应强开挖损伤区。由此，得到了巷道开挖后，掘进工作面周围煤体的系数拟合型渗透率表达式

(14)

式中，αki,βki,γki(i=1,2,3)分别为不同应力和应变条件下描述渗透率变化的系数;εip为煤进入塑性状态的应变;εir为煤进入残余状态的应变。

图4 煤的应力-应变-渗透率关系[84]Fig.4 Diagram of stress-strain-permeability for coal samples[84]

图5为掘进工作面推进过程中周围煤体内渗透率分布演化结果，周围煤体的渗透率分布同样符合三区特征。距离巷道中心较远处，煤体渗透率处于初始状态;随着与巷道中心距离的缩减，无损区内的煤体渗透率由于初始裂隙的闭合而渗透率逐渐减小(图5 Ⅰ)，距离至塑性区的边界时，渗透率达到最小值;开挖损伤区内由于产生大量裂纹而造成渗透率持续升高(图5 Ⅱ);强开挖损伤区内煤体破裂程度较为严重，由此出现大量贯通裂隙引起渗透率急剧升高(图5 Ⅲ)。对比不同开挖时间下的渗透率分布结果可知，随着时间的推移，煤体内渗透率的变化幅度呈上升趋势。

图5 掘进工作面周围煤体渗透率分布[84]Fig.5 Permeability distribution of the coal around the roadway[84]

开采扰动下工作面前方煤体或底板下煤岩体的渗透率均呈现分区段变化的特征，上述3个渗透率模型分别从断裂力学、细观损伤力学和连续损伤力学的角度考虑了采动过程中煤体损伤破裂造成的渗透率分区特征。在此之中，有效应力型的模型运用断裂力学理论，确定了微裂纹的开裂判据，通过损伤理论构建了微裂纹和损伤变量的关系，建立的模型较好的体现了采动应力和损伤破裂特征，模型的理论推导过程较为严密，具有较好的基础理论支撑。几何参数型的模型加入了煤体参数的微观统计分布，更加贴近煤体内部孔裂隙结构导致的非均匀性特征，同时模型采用拉伸和剪切破坏准则作为单元失效判据更为全面，模型从细观损伤本构出发考虑损伤破裂的因素对煤体变形和渗透的影响，适用范围较广，但计算中将煤体划分为细观单元后，求解计算的工作量较大，所以较为适合进行数值计算。系数拟合型的模型从形式上来看较为简洁，可显著区分工作面前方煤体的损伤分布特征，可以较好的反映各损伤区煤体渗透率的变化规律，此外该模型也易于进行理论计算，尤其是便于解析计算。

4 结论与展望

(1)相比天然气和煤层气的开采过程，采动应力环境下的煤层在受力方式和结构形式上更为复杂。工作面推进过程中采动应力的3个应力主方向出现不同的变化趋势，采动应力可视为外部边界条件，而煤层由于沉积作用出现的层理和割理结构属于内部结构条件。煤炭开采过程中，其内部结构和外部边界条件的主方向并非一致，所以需要考虑外部边界条件和内部渗透率主方向的夹角问题，统一两者之间的关系。与此同时，割理结构的各向异性也造成煤体物理力学参数存在异性特征。因此，需要结合地质条件和各向异性的物理力学参数建立一个考虑采动应力、煤层产状和渗透率方向的模型，以此确定煤层瓦斯渗透的最优化方向，指导煤层增透和瓦斯抽采工作。

(2)采动卸压是一个打破原始地应力后应力持续重分布引起煤层损伤破裂的增透过程。现有的煤体渗透率模型已经开始逐渐关注到这个问题，然而在已有成果中完全基于损伤理论推导出的损伤渗透率模型还较少，部分经验模型的实用性较好，但推导过程存在较多前提假设。因此，损伤渗透率模型的理论基础仍需进一步夯实。除此之外，结合断裂力学理论对煤体内裂隙演化过程进行定量描述，探究裂隙连通度对煤体渗透性的影响，确定采场尺度下不同裂隙网络空间分布与煤层渗透性的关系[85-88]，建立一个多理论交叉、宏细观结合的渗透率模型是十分必要的。

(3)煤体吸附应变的研究中，大量的实验和理论成果都是基于等温条件获得，而采面前方煤体内的气体压力和温度由表及内却是不断升高的。此外，在瓦斯吸附和解吸过程中，煤体应变严格来讲并不具有完全可逆的时空演化特征[89-91]。因此，加强非等温条件下瓦斯吸附特性的研究，掌握瓦斯解吸和吸附过程中煤体应变的差异性对采动影响下煤层渗透率演化过程的精确描述大有裨益。

总体来说，多因素影响下的采动煤体渗透率模型的理论体系并未完全建立，能够较为系统和完善体现上述特征的采动煤体渗透率模型也还鲜有出现，但煤体渗透率模型在各向异性特征、损伤破裂特征和吸附解吸特征的方面均已取得大量研究成果，以这些成果为基础的采动煤体渗透率模型将会迅速发展。同时，以上各方面的研究深度及广度仍需进一步加强。除此之外，采动应力下煤体渗透率模型的研究还有一些问题需要面对。如今，煤炭开采深度的迅速增大给煤体渗透率的研究带来巨大挑战。一方面，深部煤体的赋存环境具有“三高”特征，即高温、高地应力和高岩溶水压[92-94]。深部开采是多场多相环境共同影响的工程活动，采掘和增透过程中“三高”因素的变化相互耦合，任一因素的变化都将严重影响煤层的透气性效果。因此，深部采动煤体渗透率模型需要解决多场多相耦合的问题。另一方面，深部煤体开采过程中具有显著的强扰动和强时效特征[95-96]。强扰动是指采掘活动对深部高应力区岩体形成的强烈卸荷效应;强时效是指深部环境下采动岩体具有与采掘活动无关的、明显的流变效应。显而易见，强扰动和强时效特征将会强烈影响深部煤体的渗透特性。所以，为了全面反映深部煤体渗透率的控制因素，深部采动煤体渗透率模型还需要考虑强扰动和强时效特征的影响。