煤是天然的多孔介质，其内部含有大量破碎煤块与多种类型的微孔隙结构，且不同煤阶煤的微孔分布差异很大[1-3]。煤的微孔隙中存在孔喉结构，BAE Junseok指出煤中孔喉位置存在势垒，会对吸附/解吸过程中甲烷的通过性产生影响，导致煤中存在大量封闭孔，气孔口处的挥发性烃是微孔隙孔喉位置存在能垒的原因[4]。HE Lilin等研究表明煤的总孔隙度为7%～13%，而闭孔孔隙体积从总孔隙体积的13%～36%不等。闭孔孔隙的体积分数与显微组分组成没有相关性，但随着总孔隙体积的减小而增加[5]。Melnichenko等利用SANS与USAN技术测定了伊利诺斯盆地(美国)和Bowen Basin(澳大利亚)的3种煤中开孔的体积分数研究表明，在这些煤中的开孔的体积分数在大孔区域中的90%到中孔区域的30%之间变化，且不同煤的测试结果差异很大[6]。
孔隙的通过性对于甲烷的吸附解吸行为具有重要影响。钟玲文[7]研究表明，在常规实验条件下进行恒温解吸时，有一部分甲烷残留在煤体内部无法解吸出来;马东民等[8]指出煤层气等温解吸过程滞后于吸附过程，并引入匮乏压力下的残余吸附量对朗格缪尔方程修正。LI等研究表明，在相同解吸条件下，低阶煤残余甲烷含量小于高阶煤[9]。蔺亚兵[10]研究表明温度升高时，匮乏压力下的残余吸附量呈减小趋势，表明随着温度的增高煤层气理论解吸率和理论采收率增大。BAE Junseok指出只有通过同时测定吸附和解吸等温线，分析解吸曲线滞后特征，才能正确地估计CH4开采效率[4]。
微孔隙是煤中甲烷的主要赋存场所[11-13]，煤层气开采过程，是从微孔隙解吸，通过微孔隙喉道并进入煤中裂隙等自由空间的过程[14-15]。然而，迄今为止，甲烷在不同直径的煤的微孔隙喉道通过性及其对甲烷解吸的影响机理研究鲜见报道。因此笔者对此开展深入研究，研究结果对于煤层气开采效率评估与原位改性开采煤层气工艺的推广具有理论指导意义。
1 理论分析
1.1 甲烷在孔喉附近的势能变化特征
根据IUPAC定义孔的分类，孔径小于2 nm的孔为微孔，小于0.7 nm的孔为超微孔，甲烷主要赋存在煤的微孔中[16]。微孔的孔喉是煤指煤的微孔隙中孔径最窄的位置，甲烷在微孔隙中运移时，孔喉的煤分子与甲烷分子之间存在范德华作用力[17]，对甲烷分子运移产生明显的阻滞作用。Leonard-Jones势函数是描述分子间作用力的重要手段。煤分子与甲烷分子的Leonard-Jones势函数Ep可表示为
(1)
其中，εa为煤与甲烷分子间的吸附能，参考LIU等[18]计算出典型碳基大分子对甲烷分子的吸附能，εa 在4～15 kJ/mol;r为煤与甲烷分子间距离;r0为甲烷直径。对于二维孔喉模型(图1)，孔喉尖端的Van der Waals 径向力为
(2)
式中，R为孔喉半径;D为甲烷分子到孔喉平面的距离。
图1 甲烷分子在不同孔径孔喉附近近场势能函数
Fig.1 Adsorption potential field for different ratio of micropore mouth radius to methane molecular diameter
根据对称性可算出合力为Ftotal:
(3)
对位移积分，可以得到对应的势函数为Ep,total:
(4)
其中，x为甲烷分子相对于孔喉坐标;取εa=10 kJ/mol，可以看出:
(1)当r0≪R或r0≪|x|，即孔喉半径远大于甲烷直径、或甲烷分子与孔喉平面距离或远大于甲烷直径时，甲烷分子的势能Ep,total为0。
(2)当r0～R～|x|，即孔喉半径与甲烷直径接近时，不同孔喉半径与甲烷分子直径比的孔喉附近吸附势场如图1所示，当孔喉直径较大时，以R/r0=1.0,1.2，1.4，1.6为例，甲烷分子势能逐渐降低，并在孔喉平面出现势能最小值;随着孔喉直径进一步减小，以R/r0=0.96，0.86为例，甲烷分子势能先减小后增大，势能最小值出现在孔喉平面两侧，孔喉平面位置势能逐渐增大，甚至大于0。
1.2 甲烷分子在孔喉的通过特征分析
甲烷分子在孔喉的通过特征可基于两能态模型分析，即，设距离孔喉无穷远处的甲烷分子的势能为0;甲烷分子在通过孔喉时的最大势阱/势垒深度为Ep,max。不考虑甲烷分子之间的相互作用，单个甲烷分子通过孔喉所克服的孔喉势阱/势垒ε(R)可表示为
ε(R)= Ep,max
(5)
甲烷在孔喉位置的运移方式为碰撞-穿过，即当吸附解吸甲烷时，大量甲烷分子在孔喉平面附近进行无规则运动，并与孔喉平面发生碰撞。因此，当甲烷分子的平动动能在垂直于孔喉平面方向上的分量E大于孔喉对甲烷分子的垒时，甲烷分子可能穿越孔喉，反之，则无法通过。即:每个与孔喉平面发生碰撞甲烷分子都可能有两种状态:通过孔喉、未通过孔喉。
(6)
在R/ r0取不同值条件下对式(4)对求解最大值与最小值，可获得不同孔喉半径与甲烷直径比的势阱/势垒。如图2所示，当R/r0<0.89时，甲烷在通过孔喉时同时存在势阱与势垒，且随着孔喉减小，势阱深度接近-20 kJ/mol，而势垒迅速增大，显著大于甲烷分子在室温(278 K)条件下的平均动能(3.7 kJ/mol)，因此，孔喉对甲烷通过的阻滞作用非常大;当R/r0>0.89时，甲烷在通过孔喉时仅存在势阱，且随着孔喉增大势阱深度逐渐降低，孔喉对甲烷的阻滞作用逐渐减小;当R/r0>3时，甲烷通过孔喉时的势阱深度绝对值小于0.01 kJ/mol，显著小于甲烷分子在室温条件下的平均动能，此时，孔喉对甲烷的阻滞作用几乎消失。
图2 甲烷在不同孔径孔喉附近的势阱和势垒
Fig.2 Potential well /barrier of methane molecules at the micropore mouth with diameters
1.3 甲烷在微孔隙喉道通过性的影响因素
假设在吸附压力p(MPa)条件下甲烷吸附/解吸过程中的单位时间内，数量为N (mol)的近独立甲烷分子与半径为R(m)的圆形孔喉平面发生碰撞，由于气体分子在单位时间内对单位面积器壁的碰撞次数正是气体压强的微观表现，N与孔喉面积S，以及微孔隙内压力pinner与吸附压力p的压力差(Pa)成正比，即
N=k0ΔpS=k0πR2(pinner-p)
(7)
其中，k0为比例系数，根据麦克斯韦速度分布律，
其中k为玻尔兹曼常数，m为气体分子质量，T为温度。在热平衡条件下，与孔喉平面发生碰撞的甲烷能量服从玻尔兹曼分布[18]，能量大于某一定值ε的分子中的比例为
(8)
则得到单位时间甲烷分子的通过量Na:
(9)
由于吸附/解吸通常为恒压外部环境，若微孔内部外部的游离态甲烷服从理想气体方程
pinnerV=NinnerkT
(10)
其中，V为微孔隙体积;Ninner为孔隙内甲烷分子量。则单位时间甲烷分子的通过量也可表示为
(11)
其中，解吸时Na为正，吸附时Na为负，平衡时Na为0。
煤吸附解吸甲烷过程中，由于煤体微孔隙表面能的改变，煤体结构会发生变化，从而引起微孔隙孔喉半径的影响。即
(12)
其中，λ为变形量。煤在吸附解吸甲烷过程中的变形非常复杂。一方面，微孔中甲烷吸附/解吸时会引起微孔内、外压差的变化，导致孔喉直径发生变化;另一方面，煤吸附/解吸甲烷时，煤体骨架的膨胀与挤压变形也会引起煤微孔隙孔喉半径的改变[19-20];二者耦合导致甲烷气体在微孔中的通过性发生复杂变化。
可以看出，温度，压力差、孔喉直径是甲烷通过孔喉的主要影响因素。温度越高，压力差越大，孔喉直径越大，孔隙甲烷的通过性越好，反之则通过性越差。
2 含孔喉微孔隙解吸甲烷数值模拟
2.1 微孔隙喉道解吸甲烷数值模型
煤中甲烷的解吸是煤体从孔隙内部解吸，然后穿越孔喉进入自由空间的过程，吸附过程与之相反。为了对含微孔隙喉道煤体吸附/解吸甲烷进行数值模拟，现建立如下数值模型:
(1)把煤结构简化为大量体积为V，表面积为S的微孔隙腔，如图3所示，每个孔隙腔仅有1个孔喉半径为R的圆形孔喉。
图3 煤体结构简化为含有大量孔喉的微孔隙腔
Fig.3 Simplified coal structure models with large number micropores contains pore mouths
(2)假定微孔隙内表面存在数量为0.01 mol深度均匀的吸附势阱，且每个吸附势阱最多只吸附1个甲烷分子。如图3所示，吸附态服从朗格缪尔公式(式(13))。
其中,
(13)
其中,Nads为微孔隙内吸附态甲烷吸附量;a为吸附势阱数量;b为吸附速率参数;ε为微孔隙内表面吸附势阱深度;bm为比例系数。参考已有研究[7-8,10]，对于势阱深度为10 kJ/mol的煤体，假定在在298 K，8 MPa条件下孔隙内甲烷覆盖率达80%，代入式(13)可对比例系数bm参考值进行计算。
(3)此处不考虑煤体变形对孔喉直径的影响，即穿越孔隙的甲烷分子数服从式(11)。孔喉势阱(势垒)深度根据式(4)计算，其中εa取10 kJ/mol，r0取0.414 nm。
(4)孔隙内甲烷为吸附态与游离态两部分，即
(14)
参考已有研究[7-8]，对于吸附热(势阱深度)为10 kJ/mol的煤体，假定在在297 K，8 MPa条件下孔隙内吸附态甲烷含量约占甲烷总量的90%，根据式(13)和式(14)，可计算得到孔隙体积。
2.2 含喉道微孔隙吸附/解吸甲烷数值模拟方法
孔喉对甲烷运移的阻滞作用会对微孔隙吸附解吸甲烷速率产生显著的影响，为了完整揭示含喉道微孔隙的吸附/解吸甲烷规律，吸附/解吸实验时间均设定为5 000 min。依据微孔隙喉道解吸甲烷数值模型，对孔喉直径、微孔隙初始内部孔隙压力与外部的吸附压力，以及系统温度进行设定，然后分别对不同吸附/解吸时刻的单位时间内孔喉甲烷净通过量(mol)、微孔隙中游离甲烷量(mol)与吸附态甲烷量(mol)进行计算，从而获得微孔隙中甲烷含量随吸附/解吸时间的变化曲线与吸附解吸平衡的甲烷含量。数值模拟利用Matlab编程实现。
3 孔喉对微孔隙吸附/解吸甲烷的影响
为了研究微孔隙孔喉对微孔隙吸附解吸甲烷速率的影响，根据1.2部分分析结果，实验所选取的孔喉半径与甲烷分子直径比范围为1.09～1.57(即孔喉直径范围为0.9～1.3 nm)的孔喉。设定吸附压力为4 MPa，在303 K温度下进行微孔隙在孔喉阻滞作用下和无阻滞作用下(即在任何孔喉直径下通过率均设为100%)吸附/解吸甲烷试验。实验结果如图4所示。吸附/解吸过程中，微孔隙中甲烷吸附/解吸量与试验最终吸附/解吸量的差异小于0.1%时，则视为甲烷含量几乎不再发生变化，即达到吸附/解吸平衡。由于设定不同微孔隙吸附甲烷能力与吸附/解吸甲烷条件均一致，通过对比吸附/解吸的平衡时间即可获知孔喉对微孔隙吸附解吸甲烷速率的影响。
图4 不同孔径微孔隙孔喉对甲烷吸附/解吸的影响
Fig.4 Methane adsorption and desorption of micropores with different diameter of pore mouths
在试验时间内，无孔喉阻滞作用的孔隙(图中*所示)吸附与解吸速率极快，孔喉直径大于1.15 nm的孔隙在700 min内均达到吸附平衡，在1 800 min内达到解吸平衡;孔喉阻滞作用下的微孔隙吸附/解吸速率显著慢于无孔喉阻滞的微孔隙，孔喉直径为1.30,1.25,1.20与1.15 nm的微孔隙吸附甲烷平衡时间分别为1 300，1 800，2 700，4 500 min，小于1.15 nm的孔喉解吸甲烷平衡时间均大于5 000 min，在解吸过程中，微孔隙中甲烷的解吸速率随孔喉的减小而减小，且不同孔隙在5 000 min解吸时间内均无法完全解吸，孔喉直径越小，解吸率越低。值得注意的是，孔喉直径减小至0.9 nm(R/r0=1.09)时，微孔隙吸附与解吸甲烷的平衡时间已趋于无限长;而该孔喉直径下仅有势阱作用，无势垒作用;这表明微孔隙吸附/解吸甲烷速率缓慢是由孔喉势阱的阻滞作用引起的。
大量实验证明，煤是含有大量不同尺度微孔隙的多孔介质，煤吸附/解吸甲烷的过程不能完全可逆，即在同一吸附压力下，解吸过程中的煤中含气量始终大于吸附过程，且煤中甲烷不能够完全解吸，这种现象称为煤中甲烷解吸滞后现象[9]，该现象可导致煤层气开采过程中匮乏压力下的产率不足，导致煤层气开采效率降低。
为了阐明微孔隙孔喉与煤中解吸甲烷滞后的关系，现依据含孔喉微孔隙吸附/解吸甲烷数值模型分别构建吸附位数量为0.01 mol的4种不同微孔隙孔喉分布类型的煤样模型(表1)。根据已有研究[21-22]，不同煤阶煤体在8 MPa吸附压力下均接近饱和吸附，因此本文的甲烷吸附/解吸数值模拟主要针对小于8 MPa的压力范围进行，旨在完整揭示煤解吸甲烷滞后的机理。获得吸附解吸曲线如图5所示。可以看出，不同孔喉分布特征的煤样吸附解吸曲线均不能重合，在相同吸附压力下解吸过程含气量大于吸附过程含气量，即存在解吸滞后现象，且解吸滞后程度与孔喉直径分布特征关系密切，煤样平均孔喉直径越小，吸附解吸曲线之间区域的面积越多，煤样解吸甲烷滞后特征越明显。
表1 煤样模型
Table 1 Coal sample model
煤样编号孔喉孔径占比/%1.00 nm1.05 nm1.10 nm1.15 nm平均孔喉直径/nmNo.1402020201.06No.2204020201.07No.3202040201.08No.4202020401.09
图5 不同孔喉分布煤样模型吸附解吸甲烷曲线
Fig.5 Adsorption and desorption curves of different coal models
综上，孔喉势阱对甲烷分子运移的阻滞作用，是引起煤层甲烷吸附/解吸速率下降与甲烷解吸滞后现象，导致煤层气开采期限内解吸率低下的重要原因之一;微孔隙孔喉越小，其影响越明显。因此，通过SANS，USAN以及氮吸附等方法对储层微孔隙孔喉分布特征进行测定，是预测煤层气开采效率和产率的关键步骤。此外，通过致裂手段改变煤中孔喉的分布特征，改善甲烷分子在煤中孔喉的通过性，对于煤层气高效率开发意义重大。
4 结 论
(1)R/r0<0.89时，孔喉平面位置对甲烷分子存在势垒，孔喉平面两侧存在势阱;当R/r0>0.89时，孔喉平面仅存在势阱，势阱深度随着孔喉增大逐渐降低。
(2)温度，压力差、孔喉直径是影响甲烷通过孔喉的主要因素。温度越高，压力差越大，孔喉直径越大，孔隙甲烷的通过性越好，反之则通过性越差。
(3)孔喉势阱对甲烷分子运移的阻滞作用，是引起煤层甲烷吸附/解吸速率下降与甲烷解吸滞后现象，导致煤层气开采期限内解吸率低下的重要原因之一;微孔隙孔喉越小，其影响越明显。
[1] CHENG Yuanping,JIANG Haina,ZHANG Xiaolei,et al.Effects of coal rank on physicochemical properties of coal and on methane adsorption[J].International Journal of Coal Science & Technology,2017,4(2):129-146.
[2] LIU Z,ZHANG Z,CHOI S K,et al.Surface properties and pore structure of anthracite,bituminous coal and lignite[J].Energies,2018,11(6):1502.
[3] 蔺亚兵,贾雪梅,马东民,等.基于液氮吸附法对煤的孔隙特征研究与应用[J].煤炭科学技术,2016,44(3):135-140.
LIN Yabing,JIA Xuemei,MA Dongmin,et al.Study and application of coal pore features based on liquid nitrogen adsorption method[J].Coal Science & Technology,2016,44(3):135-140.
[4] BAE J S,BHATIA S K,RUDOLPH V,et al.Pore accessibility of methane and carbon dioxide in coals[J].Energy & Fuels,2009,23(6):3319-3327.
[5] HE L,MELNICHENKO Y B,MASTALERZ M,et al.Pore accessibility by methane and carbon dioxide in coal as determined by neutron scattering[J].Energy Fuels,2012,26(3):1975-1983.
[6] MELNICHENKO Y B,HE L,SAKUROVS R,et al.Accessibility of pores in coal to methane and carbon dioxide[J].Fuel,2012,91(1):200-208.
[7] 钟玲文.煤的吸附性能及其影响因素[J].地球科学,2004,29(3):238-332.
ZHONG Lingwen.Adsorption properties of coal and its influencing factors[J].Earth Science,2004,29(3):238-332.
[8] 马东民,张遂安,蔺亚兵.煤的等温吸附-解吸实验及其精确拟合[J].煤炭学报,2011,36(3):477-480.
MA Dongmin,ZHANG Suian,LIN Yabing.Isothermal adsorption desorption experiment of coal and its accurate fitting[J].Chinese Journal of Coal Science,2011,36(3):477-480.
[9] LI Z,LIN B,GAO Y,et al.Fractal analysis of pore characteristics and their impacts on methane adsorption of coals from Northern China[J].International Journal of Oil Gas & coal Technology,2015,10(3):306.
[10] 蔺亚兵.煤层气解吸滞后效应研究[D].西安:西安科技大学,2012.
LIN Yabing.Study on desorption hysteresis effect of coalbed methane[D].Xi’an:Xi’an University of Science and Technology,2012.
[11] KANG J,ZHOU F,XIA T,et al.Numerical modeling and experimental validation of anomalous time and space subdiffusion for gas transport in porous coal matrix[J].International Journal of Heat & Mass Transfer,2016,100:747-757.
[12] PANDEY R,HARPALANI S,FENG R,et al.Changes in gas storage and transport properties of coal as a result of enhanced microbial methane generation[J].Fuel,2016,179:114-123.
[13] CUI X,BUSTIN R M.Volumetric strain associated with methane desorption and its impact on coalbed gas production from deep coal seams[J].AAPG Bulletin,2005,89(9):1181-1202.
[14] YANG X,HUANG F,DAI W,et al.Thermodynamic performance analysis for flashing-binary cycle using different working fluids to low temperature steam[J].Proceedings of the CSEE,2014,34(20):3257-3265.
[15] QIN L,ZHAI C,LIU S,et al.Failure mechanism of coal after cryogenic freezing with cyclic liquid nitrogen and its influences on coalbed methane exploitation[J].Energy & Fuels,2016,30(10):8567-8578.
[16] YANG H P,MING L I.Study on the adsorption force field of methane on coal surface molecules[J].Journal of Xi’an University of Science & Technology,2015,35(1):96-99.
[17] 聂百胜,段三明.煤吸附瓦斯的本质[J].太原理工大学学报,1998(4):417-421.
NIE Baisheng,DUAN Sanming.Essence of coal adsorbed gas[J].Journal of Taiyuan University of Technology,1998(4):417-421.
[18] LIU Z,FENG Z,ZHANG Q,et al.Heat and deformation effects of coal during adsorption and desorption of carbon dioxide[J].Journal of Natural Gas Science & Engineering,2015,25:242-252.
[19] FENG Z,ZHOU D,ZHAO Y,et al.Study on microstructural changes of coal after methane adsorption[J].Journal of Natural Gas Science & Engineering,2016,30:28-37.
[20] 
R,THOMAS H R.Effect of sorption induced swelling on gas transport in coal[A].GeoShanghai International Conference[C].Springer,Singapore,2018:386-394.
[21] ZHANG Qingling,CUI Yongjun,CAO Lige.Effect of pressure on adsorption performance of coals with different degrees of metamorphism[J].Natural Gas Industry,2004,24(1):98-100.
[22] 马东民,马薇,蔺亚兵.煤层气解吸滞后特征分析[J].煤炭学报,2012,37(11):1885-1889.
MA Dongmin,MA Wei,LIN Yabing.Analysis of desorption hysteresis characteristics of coalbed methane[J].Journal of Coal Science,2012,37(11):1885-1889.
