沁水盆地南部无烟煤大分子结构模型及其含甲烷力学性质

张 彬，曾凡桂，王德璋，康官先，张晓雨，康天合

(1.太原理工大学 原位改性采矿教育部重点实验室，山西 太原 030024; 2.太原理工大学 矿业工程学院，山西 太原 030024; 3.山西晋城无烟煤矿业集团有限责任公司，山西 晋城 048006; 4.太原理工大学 煤与煤系气地质山西省重点实验室，山西 太原 030024; 5.太原理工大学 安全与应急管理学院，山西 太原 030024)

沁水盆地南部的煤层气开发已实现了小型商业化，但仍存在产量低和衰减快的问题[1-2]，往往需要压裂改造提高抽采率[3]，而煤层本身的力学性质是决定压裂效果的关键因素[4]。陈立超和王生维[5]利用测井数据反演了郑庄区块煤岩的弹性模量、剪切模量、体积模量和泊松比等参数。冯晴等[6]实验得出弹性模量和泊松比随埋深的变化规律。实际上，沁水盆地无烟煤储层是普遍含气的[7]，含甲烷条件下的煤体力学性质更符合生产实际。辛程鹏等[8]对新景矿含瓦斯突出煤进行了常规三轴和分段变速加载力学试验研究，结果表明:随瓦斯压力降低，煤体强度和弹性模量均增大;煤体在分段变速加载路径下的强度普遍增大，峰值轴向应变、峰值环向应变绝对值和峰值体积应变绝对值也普遍增大，失稳破坏瞬间应力跌落和能量释放更加剧烈。孟筠青等[9]采用分子模拟的方法对赵庄煤大分子模型进行力学性质研究，并分析了径向分布曲线和键角分布曲线，发现煤大分子中的环状结构在抵抗破坏过程中起到的重要作用，提出了研究纳米尺度下煤的力学性质的新思路，是揭示气体与固体相互作用机制的有力手段[10]。笔者以沁水盆地15号煤为研究对象，对其结构进行表征并构建其大分子结构模型。在模型的基础上，采用巨正则系综蒙特卡洛(GCMC)对甲烷在无烟煤中的吸附特性进行研究，获得吸附量、吸附热和吸附构型;采用分子动力学的方法对吸附甲烷后的构型进行力学性质计算，从分子角度揭示含甲烷无烟煤的力学性质，为沁水盆地煤层气抽采提供一些微观的理论参数。

1 研究方法和研究方案

分子模拟研究含甲烷无烟煤力学性质的方法为:① 无烟煤化学结构表征;② 基于表征数据，建立无烟煤大分子结构模型;③ 模拟无烟煤大分子结构模型吸附甲烷分子;④ 计算分析含甲烷无烟煤构型的力学性质。

1.1 无烟煤化学结构表征

目前，常采用工业分析、元素分析、核磁共振碳谱(13C-NMR)和X射线光电子能谱(XPS)等手段来表征煤体结构，如马汝嘉等[11]表征了陕西凤县高煤级煤的化学结构。笔者也采用相似的方法对沁水盆地南部无烟煤进行表征。

表1为无烟煤工业分析和元素分析的测试结果。测试依据为《煤的工业分析方法》(GB/T212—2008)和《煤的元素分析方法》(GB/T31391—2015)进行无烟煤样品的工业分析(以空气干燥基为标准)和元素分析(以干燥无灰基为标准)。

表1 无烟煤工业分析和元素分析Table 1 Proximate and ultimate analysis of anthracite sample %

图1 无烟煤样的13C-NMR频谱分峰拟合Fig.1 13C-NMR spectral peak fitting of anthracite sample

表2 13C-NMR 化学位移归属Table 2 13C-NMR chemical shift attribution

表3 无烟煤13C-NMR结构参数计算结果Table 3 Calculation results of 13C-NMR structural parameters of anthracite %

图2分别给出了无烟煤XPS总能谱以及C1s,N1s,S2p能谱范围的分峰拟合曲线。测试仪器为ESCALA B250型的X射线光电子能谱仪。利用Avantage软件分别对C,N，S元素的能谱进行分峰拟合确定O,N，S元素不同存在形式的相对含量。表4给出了O,N，S元素不同化合态的分布情况。

1.2 无烟煤大分子结构模型的构建

根据上述表征结果，得出沁水盆地15号无烟煤的分子式为C228H77N3O13S3，共由324个原子组成，其大分子模型的结构参数和芳香结构见表5，6。构建的沁水盆地无烟煤结构平面模型，如图3所示。将5个无烟煤的结构平面模型沿三维周期性组合，构建无烟煤大分子结构模型，如图4所示。无烟煤大分子结构模型的三维尺寸为2.74 nm ×2.74 nm ×2.74 nm。

图2 无烟煤XPS总能谱以及C1s,N1s,S2pFig.2 XPS spectrum and fitting curves of C1s,N1s,S2p of anthracite

表4 无烟煤杂原子分布特征Table 4 Distribution characteristics of hybrid atoms in anthracite

表5 无烟煤大分子模型的结构参数Table 5 Parameter of molecular model unit of anthracite

表6 无烟煤大分子模型的芳香结构
Table 6 Aromatic structure of macromolecular model of anthracite

图3 沁水盆地15号无烟煤的结构平面模型Fig.3 Structural plane model of No.15 Anthracite in Qinshui Basin

图4 沁水盆地无烟煤的大分子结构模型Fig.4 Macromolecular structure model of No.15 anthracite in Qinshui Basin

1.3 无烟煤吸附甲烷模拟方法和方案

目前，许多研究表明蒙特卡洛(GCMC)分子模拟方法是研究气体在多孔介质中性质的有效手段[14,15]。因此，笔者通过Materials Studio软件中的Sorption模块，采用蒙特卡洛(GCMC)法研究甲烷在无烟煤中的吸附特性[16]。模拟力场为COMPASS，方法为Metropolis，Quality为Customized，电子势和范德华势分别采用Ewald和Atom based方法进行统计处理，模拟1×107步，其中前5×106步用于使体系达到平衡，后5×106步用于计算吸附量和吸附热等热力学参数。

考虑目前煤层气资源量以埋深2 000 m为临界值，其原岩应力近似于30 MPa。因此，模拟计算时，设置压力为1,5,10,20,30 MPa，设置温度为50 ℃。在Materials Studio软件中需要输入逸度而非压力，通过Peng-Robinsion方程实现压力和逸度的转换[17]。

1.4 含甲烷无烟煤力学性质模拟方法

对固体材料而言，通常利用杨氏模量、体积模量、剪切模量和泊松比等参数来表征其力学性质。为更好的研究无烟煤的含气力学性质，基于上述蒙特卡洛法得到1,5,10,20,30 MPa下甲烷在无烟煤中的吸附构型，采用分子动力学的方法，在Materials Studio软件中Forcite模块的Mechanical Properties程序中以恒定应力的方式对构型进行单轴拉伸和纯剪切变形计算，其中最大应变为0.003，压强为1个大气压(0.1 MPa)，迭代次数为10 000。在体系始终保持稳定的前提下，材料的力学参数[18]可表示为

(1)

(2)

G=μ

(3)

(4)

式中，Y为弹性模量，GPa;λ和μ分别为Lame常数，GPa;K为体积模量，GPa;G为剪切模量，GPa;ν为泊松比。

2 甲烷在无烟煤中的吸附特性

2.1 模型的正确性验证

为验证模型的正确性，普遍将模拟数据和实验数据进行对比。图5为无烟煤平面模型谱和固体核磁谱的对比。可以看出，无烟煤平面模型的核磁谱与上述实验的固体核磁谱具有较好的一致性，对模型进行反演计算得出无烟煤的元素含量，与之前实验煤样进行元素分析时，含量的结果非常接近，与理论测定无烟煤的元素含量范围[19]相符，详见表7。采用分子动力学模拟得到的无烟煤模型密度为1.526 g/cm3，符合前人通过自然伽马测井测定的密度(1.39～1.76 g/cm3)[20]，也满足前人理论计算得到的无烟煤密度(1.3～1.7 g/cm3)[21]。根据理想气体定律，用TALU和MYERS[22]的方法模拟计算了无烟煤的孔体积:

图5 无烟煤平面模型谱和固体核磁谱的对比Fig.5 Comparison of plane model spectrum and solid state nuclear magnetic resonance spectrum of anthracite

其中，R为气体常数;Nm为单位质量吸附剂(Mm)的吸附探针氦气的分子数;T为温度;P为压力。模拟得到的无烟煤大分子结构模型的孔体积为0.018 6 cm3/g，与戚灵灵等[23]采用液氮吸附测得寺河矿煤的孔体积(0.016 4 cm3/g)相差不大。因此，通过上述分析说明本文所建立的无烟煤平面模型和大分子结构模型是合理的，能够保证后续计算的正确性。

表7 模拟计算、实验测定和理论范围的无烟煤元素质量分数、密度和孔体积Table 7 Element content,density and pore volume of anthracite in simulated calculation,experimental measurement and theoretical range

2.2 吸附量和吸附构型

图6为甲烷在无烟煤中的吸附量随压力的变化曲线。可以看出，甲烷在无烟煤中吸附量呈先增大后趋于稳定的变化规律，符合Langmuir模型[24]。压力PL为1,5,10,20和30 MPa时，甲烷在无烟煤中的吸附量VL分别为9.64,15.64,17.68,19.67,21.03个/晶胞，呈先增大后趋于稳定的变化规律，符合Langmuir模型。通过Langmuir模型对数据进行拟合，得出甲烷在无烟煤中的饱和吸附量为22.4个/晶胞，Langmuir压力为1.12 MPa。图7为20 MPa时甲烷在无烟煤中的吸附构型，其他压力条件也具有相似的构型，但是吸附甲烷数量不同。

图6 甲烷在无烟煤中的吸附量随压力的变化曲线Fig.6 Curve of methane adsorption in coal with pressure

图7 20 MPa时甲烷在无烟煤中的吸附构型Fig.7 Adsorption configuration of methane in anthracite at 20 MPa

2.3 吸附位

径向分布函数(RDFs)是描述原子密度与原子间距离的函数，可以用来反映煤等多孔介质材料中的不同原子与甲烷等气体分子的相互作用强度，判断其吸附位[25]表达式为

(5)

式中，gij(r)为在分子调围距离r处，分子j的局部密度;dN为从r到r+dr的分子j的数量;ρj为分子j的密度。

图8为压力20 MPa时，无烟煤模型中的碳、氢、氧、氮和硫原子与甲烷分子的径向分布函数曲线，可以看出第1个RDFs波峰出现在距离为0.11 nm处，为CH4—C的波峰，峰值为5.78。其次是在距离为0.238 nm和0.43 nm处分别由CH4—H和CH4—N产生的波峰，CH4—O和CH4—S引起的波峰距离较远，波峰相对较低。因此，芳香碳、吡啶型氮和吡咯型氮以及羧基是甲烷分子主要吸附位。甲烷分子更容易吸附在由芳香环和脂肪结构组成的孔隙中，且更容易被孔隙中靠近芳香结构的一侧吸附。

图8 无烟煤与甲烷分子的径向分布函数曲线Fig.8 Radial distribution function of anthracite and methane molecule

2.4 等量吸附热

在甲烷吸附过程中，等量吸附热能够反映无烟煤与甲烷间能量释放的变化关系，可通过Clausius-Clapeyron方程计算求得[26]。图9为甲烷在无烟煤中的等量吸附热随吸附压力的变化规律。可以看出，甲烷在无烟煤中的等温吸附热qst随吸附压力P的升高呈对数规律下降，拟合得到的对数表达式为qst=-0.71× lnP+25.05，表明甲烷在无烟煤中吸附过程中释放的热量随着压力的升高逐渐减小，并最终达到平衡，这是因为甲烷在低压力时率先占据无烟煤表面的高能吸附位。在吸附压力范围内，甲烷在煤中的吸附热为22.65～25.00 kJ/mol，远小于42 kJ/mol，属于物理吸附。

图9 甲烷在煤中的等量吸附热随压力的变化曲线Fig.9 Curve of isosteric heat of methane adsorption in coal with pressure

3 无烟煤的含气力学性质

3.1 力学性质

表8为不同含气量无烟煤的力学参数。图10为无烟煤的力学参数随吸附量Q的变化规律。由图10可以看出，体积模量K、弹性模量E和剪切模量G具有一致的变化规律，均随着甲烷吸附量的增大呈对数规律降低，而泊松比ν随吸附量的增大呈线性规律增大，拟合方程分别为

K=ln(881.04-36.65Q)

(6)

E=ln(234.15-8.07Q)

(7)

G=ln(7.71-0.14Q)

(8)

ν=ln(0.31+0.002Q)

(9)

王晋等[27]研究了不同含气饱和度下煤岩力学性质的变化趋势，对其数据进行拟合(图11)，发现煤的弹性模量与含气饱和度呈对数衰减规律，与本文的研究结论相一致。当甲烷含气量为21.03个/晶胞时，无烟煤的体积模量、弹性模量和剪切模量分别为4.12,4.29和1.45 GPa，同比于不含甲烷的无烟煤，体积模量、杨氏模量和剪切模量分别降低了38.5%,24.4%和27.1%。这表明吸附甲烷能够显著降低无烟煤的力学强度，且吸附量越大影响程度越高。

图10 无烟煤的力学参数随吸附量的变化曲线Fig.10 Curves of the mechanical parameters of coal with the adsorption capacity

图11 煤的弹性模量随含气饱和度的变化曲线Fig.11 Variation curve of elastic modulus of coal with gas saturation

3.2 结果讨论与机理分析

高保彬等[28]采用自主研发的含瓦斯煤岩受载变形破坏试验装置，通过实验室试验研究了不同瓦斯压力条件下无烟煤的力学性质，设置瓦斯压力为0～2 MPa，发现煤的力学强度随瓦斯压力的增大而降低，与本文的结论相一致。但由于其所设置压力范围偏小，以此结果进行无烟煤力学性质评价有一定的局限性。本试验延伸了瓦斯压力范围(20 MPa)，定量分析了无烟煤力学性质和吸附量之间的关系，得出了无烟煤的体积模量、杨氏模量和剪切模量随着吸附量的变化规律。

已有研究表明，多孔介质材料吸附甲烷后会发生膨胀变形，其计算公式[16]为

a=(Vi-V0)/V0

(10)

其中，a为膨胀率;Vi，V0分别为煤的原始体积和吸附甲烷后的体积。图12为无烟煤的体积和膨胀率随含气量的变化规律。由图12可知，无烟煤的体积和膨胀率随吸附量的增加呈指数规律增大。与未吸附甲烷的无烟煤相比，当甲烷含气量为21.03个/晶胞时，无烟煤的膨胀率达到了9.012%。根据周动等[29]的研究，发现无烟煤的膨胀变形量随煤结构尺度的减小而显著增大。当煤结构尺度减小时，膨胀变形量的降幅增大，这是由于在较小统计尺度下煤结构的非均质性强，结构间力学强度与膨胀变形能力差异较大，膨胀变形量均较大;而在较大统计尺度下煤结构的非均质性弱，结构间力学强度与膨胀变形能力差异较小，膨胀较小。煤体的膨胀变形主要受到吸附压力与结构尺度的双重影响。ZHANG等[30]采用分子模拟的方法研究了含水烟煤吸附甲烷后的膨胀变形，发现当含水率为3.0%和压力为1 MPa时，烟煤的膨胀率达到3.1%。本文的研究是不考虑含水状态下的含甲烷力学性质，甲烷的吸附位没有被水占据，能够吸附更多的甲烷，产生更大的膨胀变形。因此，本文在纳米尺度条件下，无烟煤吸附甲烷的膨胀率是可以接受的。

表9为不同含气量无烟煤模型的能量对比。可以看出，随着含气量的增大，煤体会发生吸附膨胀，导致无烟煤模型的总能量降低，说明无烟煤基质之间的相互作用力降低，进而导致无烟煤的强度降低，抵抗变形的能力减弱。表9结果表明，虽然成键能与非成键能相比降低幅度小，但键能、键角能、扭转能、反转能等均有不同程度的降低，这说明甲烷的吸附将导致键能的变化，使其键强降低。分子间作用力是物质稳定的重要物理因素，表中数据表明范德华能在吸附甲烷的过程中能量降低值最大，对无烟煤总能量的降低起主导作用，表明范德华能是保持煤体结构和力学性质稳定的主导因素。当含气量分别为9.64，15.64，17.68，19.67和21.03个/晶胞时，无烟煤的范德华能分别降低了25.02，27.97，28.99，30.86和32.60 kJ/mol。因此，含甲烷无烟煤体积模量、弹性模量、剪切模量等降低、泊松比增大的内在原因是吸附甲烷会降低煤分子的键能和分子间相互作用力。这种内在关系值得进一步的理论与实验研究。

图12 无烟煤的体积和膨胀率随吸附量的变化曲线Fig.12 Variation of volume and expansion rate of anthracite with gas content

表9 不同含气量无烟煤模型的能量对比Table 9 Energy comparison of anthracite models with different gas contents kJ/mol

4 结 论

(1)无烟煤样大分子结构中主要由芳烃类碳组成，平均缩聚芳香结构由3～4个苯环的链型连接或环型连接。

(2)甲烷在无烟煤中吸附量呈先增大后趋于稳定的变化规律，符合Langmuir模型，甲烷在无烟煤中的饱和吸附量为22.4个/晶胞，Langmuir压力为1.12 MPa。甲烷分子更容易吸附在由芳香环和脂肪结构组成的孔隙中，且更容易被孔隙中靠近芳香结构的一侧吸附。

(3)无烟煤的体积模量、弹性模量和剪切模量均随着甲烷吸附量的增大呈对数规律降低，而泊松比随吸附量的增大呈线性增大。

(4)随着含气量的增大，煤体会发生吸附膨胀，无烟煤模型的总能量降低，表明无烟煤基质之间的相互作用力降低，进而导致无烟煤的强度降低，抵抗变形的能力减弱。范德华能是保持煤体结构和力学性质稳定的主导因素。

(5)分子模拟是研究煤等多孔介质材料在含气状态下力学性质的一种有效手段。