伴随世界油气需求的持续增长与常规油气产量的不断下降，具有较大资源潜力的非常规油气逐渐成为勘探开发的热门领域，受到各个国家和石油公司的高度重视。页岩储层中气体的赋存状态主要包括:游离态、吸附态及溶解态。据估计，吸附气占页岩气地质储量的比例为20% ～ 85%[1-2]。因此，页岩的吸附特性研究对页岩气开采具有重要意义。
在页岩储层中，页岩气主要吸附于有机质(干酪根)中和无机质(黏土矿物)表面[3]。有机质对页岩气的吸附能力较强，许多研究认为页岩吸附量与总有机碳(Total Organic Content，TOC)含量呈正相关[4-6]。与此同时，页岩中的黏土矿物对甲烷的吸附量也相对较大。1991年，SCHETTLER等[7]发现页岩中黏土矿物吸附量很大，某些情况下，页岩吸附曲线与伊利石吸附曲线相近。1995年，LU等[8]测量了Devonian页岩样品和纯伊利石的吸附能力，认为虽然TOC含量是影响页岩吸附能力的首要因素，但黏土矿物的吸附能力同样不容小觑，尤其在TOC含量较低的页岩中。2009年，ROSS和BUATIN[9]在进行页岩组成及孔隙结构研究的室内试验时发现，除有机质外，页岩中黏土矿物的微孔体积也很大。此外多数研究认为页岩中除有机质和黏土矿物以外，其他矿物中基本无微孔[5,10-11]，所以其吸附量基本可忽略[12]。
页岩中的黏土矿物以蒙脱石、高岭石、伊利石和伊/蒙混层等矿物为主[13]，美国Barnett页岩气藏中伊/蒙混层与伊利石的比例为4∶1[14]，杨峰等[15]研究中国四川盆地下寒武统牛蹄塘组和下志留统龙马溪组海相页岩样品时发现，页岩样品中黏土矿物含量平均为32.2%，且主要为伊/蒙混层和伊利石。近年来，随着计算机分子模拟手段的广泛使用，已经有学者利用分子模拟手段开展了甲烷在蒙脱石、伊利石、高岭石等页岩主要矿物上吸附行为的研究[16-22]，但目前鲜有针对甲烷在伊/蒙混层上吸附的分子模拟研究。因此笔者构建了可靠的伊/蒙混层模型，通过分子模拟手段，对比不同黏土矿物对甲烷吸附能力的差异，并从微观角度研究甲烷在伊/蒙混层中吸附行为的影响因素，为页岩气的勘探开发及储量预测提供借鉴。
1 分子模型
1.1 蒙脱石模型的构建
本文所使用的蒙脱石模型根据SKIPPER等[23]的蒙脱石模型构建。蒙脱石晶胞参数为:晶胞三组棱长的长度a=0.523 nm，b=0.906 nm，c=1.25 nm，棱间交角α=γ=90°，β=99°，分子式为Na0.75(Si7.75Al0.25)(Al3.5Mg0.5)O20(OH)4。蒙脱石的分子模型如图1所示。
图1 蒙脱石分子模型
Fig.1 Molecular model of smectite
1.2 伊利石模型的构建
本文所使用的伊利石晶体模型根据DAWN等[24]提出的几何构型构建。1M-tv伊利石模型晶胞参数为a=0.519 nm，b=0.895 nm，c=0.995 nm，α=γ=90°，β=94.87°，分子式为KAl4(OH)4(Al0.5Si3.5O10)2。伊利石的分子模型如图2所示。
图2 伊利石分子模型
Fig.2 Molecular model of illite
1.3 伊/蒙混层模型的构建
伊/蒙混层(I/S)是蒙脱石向伊利石转化过程中形成的中间产物，是由蒙脱石晶层和伊利石晶层沿C轴方向组成的特殊类型的层状硅酸盐矿物。目前认为的伊/蒙混层晶体模型有3种[25]，分别是非极性2∶1层状模型、极性2∶1层状模型以及多相模型。本文采用的伊/蒙混层模型根据ALTANER[25]的非极性2∶1层状模型构建，蒙脱石片层与伊利石片层相互间隔。伊/蒙混层的分子模型如图3所示。
图3 伊/蒙混层分子模型
Fig.3 Molecular model of mixed-layer illite/smectite
1.4 分子模拟细节
所有的模拟单元均由4×2×1的超晶胞构成，模拟过程采用周期性边界条件，将两层黏土层组成狭缝状纳米孔结构，孔径定义为内表面层中氧原子中心平面之间的距离。黏土矿物建模使用CYGAN等[26]的CLAYFF力场，甲烷分子采用TraPPE模型描述[27]。
在模拟吸附之前，对矿物模型进行结构优化，使其能量达到最小。采用Forcite模块中的Dynamics功能对黏土矿物模型进行优化，以获得稳定的构型。选取正则系综(NVT)进行500 ps的动力学模拟，模拟步长为1 fs，温度设定为338.55 K，控温系统采用Nose热浴，力场选择CLAYFF力场，原子类型和电荷由力场自动分配，静电力加和及范德华力加和方法分别为Ewald和Atom based，截断半径设置为0.9 nm。每次GCMC模拟所加载的最大步数为2×107，其中前1×107用于吸附平衡，后1×107步用于平衡吸附量的数据统计。在模拟过程中，黏土片层被假定为刚性的，因为研究的重点是吸附机制，而不是吸附引起的黏土膨胀。层间离子不固定，以使层间离子可以在力场的作用下自由移动。
2 结果与讨论
2.1 分子模型验证
INOUE等(1987)[28]通过对黏土矿物的结构进行化学分析，认为自然界中基本单元厚度为1 nm的结构属于蒙脱石结构;厚度在2～4 nm的结构介于端元蒙脱石与伊利石之间，属于伊/蒙混层结构;厚度≥5 nm的结构属于伊利石结构。以此为依据构建了孔径分别为1，2和9 nm的蒙脱石、伊/蒙混层及伊利石模型。图4为甲烷在3种黏土矿物模型中的过剩吸附量，模拟的温度和压力与文献[29]的实验条件相同。模拟结果显示，3种黏土矿物对甲烷的吸附能力顺序为:蒙脱石>伊/蒙混层>伊利石，这与实验数据显示的趋势一致。分子模拟得到的甲烷在蒙脱石层、伊/蒙混层和伊利石层中最大过剩吸附量分别是0.399，0.251，0.062 mmol/g，而文献[29]中按照黏土矿物含量折算得到的过剩吸附量分别为0.397，0.225，0.045 mmol/g，证实了模型的可靠性。
图4 338.55 K甲烷在3种黏土矿物模型中的过剩
吸附量曲线
Fig.4 Excess adsorption isotherms of methane on 3 clay mineral models under 338.55 K
模拟得到的甲烷吸附量略高于实验值，主要原因可能有:① 进行模拟的黏土矿物为无水化理想模型，而实验所用的样品可能含有吸附水，吸附水的存在会造成黏土矿物对甲烷的吸附量减小;② 用于模拟的矿物模型为单一纯矿物，而实验所用的样品中可能含有影响甲烷吸附的杂质;③ 黏土矿物模拟与实际矿物有差别，如层间阳离子、原子类质同相取代位置不同等;④ 进行吸附实验前，矿物样品被粉碎成细小颗粒，颗粒之间的孔隙属于大孔级别[30]，而模型孔隙属于微孔级别，因此模拟得到的过剩吸附量大于实验结果。
2.2 等量吸附热
模拟得到的甲烷在孔径分别为1，2，4 nm的3种上黏土矿物上的平均等量吸附热如图5所示，平均等量吸附热均小于42 kJ/mol，说明甲烷在黏土矿物表面发生物理吸附[31]。
图5 甲烷在不同孔径黏土矿物上的等量吸附热
Fig.5 Equivalent adsorption heat of methane on clay minerals with different pore sizes
吸附热越大，表明吸附能力越强，3种黏土矿物对甲烷的吸附能力为:蒙脱石>伊/蒙混层>伊利石。随着压力增加，甲烷与黏土矿物片层的接触几率增加，同时甲烷分子间的相互作用力增强，导致黏土对甲烷的吸附能力增强，吸附热增加。孔径越小，等量吸附热越大，证明吸附越牢固，这一现象可以从吸附势理论的角度得到解释。孔径较小时，相对的两个孔壁距离很近，孔壁产生的吸附势重叠[32]，对吸附质分子的作用力比大孔更大，即对气体的吸附能力增强，因此等量吸附热大。
2.3 等温吸附量
页岩中黏土矿物的孔隙结构具有很强的非均质性，因此研究孔径对页岩中黏土矿物吸附能力的影响是十分必要的。分别构建孔径为1，2，4 nm的黏土矿物模型，模拟温度338.55 K、压力1～20 MPa条件下甲烷在3种模型上的吸附行为，模拟结果如图6，7所示。
图6 338.55 K甲烷在不同黏土矿物上的吸附等温线
Fig.6 Adsorption isotherms of methane on different clay minerals at 338.55 K
分子模拟得到的吸附量为总吸附量，与过剩吸附量之间的转换关系为
nex=nab-Vg ρg
(1)
式中，nex为气体过剩吸附量，mmol/g;nab为气体总吸附量，mmol/g;Vg为黏土矿物的孔隙空间，cm3/g;ρg为流体体相密度，mmol/cm3，由Peng-Robinson EOS方程[33]求得。
模拟结果表明，甲烷在3种黏土矿物上的吸附呈相同规律。甲烷总吸附量随压力升高而缓慢增加，最终趋于恒定，孔径越小，总吸附量越小，达到饱和吸附量所需的压力越低。原因是甲烷在黏土矿物上发生多层吸附，孔径空间较小时，能够进入孔隙的气体分子数量有限，吸附层数受到有限空间的限制，因此吸附量小。
随着压力升高，甲烷过剩吸附量逐渐增加，达到最大值后下降。根据式(1)，产生这一现象的原因是甲烷过剩吸附量受到总吸附量和甲烷体相密度的影响，甲烷体相密度随压力增加而缓慢增大。低压阶段，甲烷吸附相密度大于体相密度，此时过剩吸附量为正值，且过剩吸附量随着压力升高而逐渐増加，在吸附相密度与体相密度相等时对应的压力点处达到极值，此后随着压力继续升高，吸附相密度小于体相密度，导致过剩吸附量为负值。孔径越小，吸附相密度越大，因此过剩吸附量越大。
图7 338.55 K甲烷在不同黏土矿物上的过剩吸附量
Fig.7 Excess adsorption of methane on different clay minerals at 338.55 K
CH4在黏土矿物表面的吸附等温曲线符合I型吸附曲线特征，采用Langmuir等温吸附模型进行拟合，用吸附量表示的Langmuir方程为
(2)
其中，x为气体的绝对吸附量，mmol/g;am为气体的饱和吸附量，mmol/g;k为Langmuir平衡常数，MPa-1;p为压力，MPa。拟合结果见表1，可见Langmuir模型对CH4吸附数据的拟合精度较高。
2.4 甲烷在伊/蒙混层狭缝中的空间分布
为了进一步研究甲烷分子在黏土矿物狭缝中的空间分布特征及数量变化关系，绘制甲烷分子沿垂直于孔壁方向的密度分布曲线。鉴于目前鲜有针对甲烷在伊/蒙混层上吸附的分子模拟，本文以伊/蒙混层为例进行瞬时构型描述。338 K，20 MPa下，甲烷分子在不同孔径伊/蒙混层狭缝中的瞬时空间分布和密度剖面如图8和9所示。可以看出，相同条件下，甲烷气体在孔径为1 nm的狭缝中形成2个吸附层，密度峰值较高，此时甲烷仅以吸附态存在;在大于2 nm的狭缝中，随着孔径增加，吸附层数逐渐增多，吸附相密度逐渐降低，此时甲烷以吸附态和游离态共同存在。在不同孔径的孔隙中，蒙脱石一侧甲烷吸附层的密度始终高于伊利石一侧，表明伊/蒙混层矿物中蒙脱石片层的吸附能力强于伊利石片层。
表1 甲烷吸附等温线的Langmuir拟合结果
Table 1 Langmuir fitting results of methane adsorption isotherm
孔径/nm蒙脱石am/(mmol·g-1)k/MPa-1R2伊/蒙混层am/(mmol·g-1)k/MPa-1R2伊利石am/(mmol·g-1)k/MPa-1R212.6610.0760.9992.1700.091 70.9992.1570.0890.99927.7820.0300.9987.5190.029 40.9987.2360.0290.999420.2020.0180.99819.9600.017 90.99919.1570.0180.999
温度为338 K，不同压力条件下，甲烷在2 nm伊/蒙混层狭缝中的瞬时空间分布和密度剖面如图10和11所示。可以清晰地看出，黏土矿物狭缝中，层间阳离子比甲烷分子更靠近黏土片层。低压情况下，层间阳离子占据黏土片层附近的吸附位，而甲烷占据离子旁边的位置;压力逐渐升高，气体分子开始占据黏土矿物附近的空位。吸附气量随着压力增加逐渐增大，低压条件下，黏土矿物壁面的吸附作用是影响气体吸附量的主要因素;但随着压力升高，压力的作用逐渐占据主导地位，甲烷气体在高压下受到压缩，密度增加，导致气体之间的相互作用增大，相比之下吸附作用的影响十分微小。压力较低的情况下，蒙脱石和伊利石片层上吸附气的密度差异较为显著，随着压力升高，两侧密度差异逐渐减小，表明吸附作用对吸附相密度的影响随压力增加而减弱。
图8 甲烷分子在不同孔径伊/蒙混层中的瞬时构型
Fig.8 Snapshots of methane adsorption on mixed-layer illite/smectite surfaces,pore sizes
图9 甲烷在不同孔径伊/蒙混层中的密度剖面
Fig.9 Density profiles of methane adsorption on mixed-layer illite/smectite,pore sizes
图10 不同压力下，甲烷分子在孔径2 nm伊/蒙混层中的
瞬时构型
Fig.10 Snapshots of methane molecule adsorption on 2 nm mixed-layer illite/smectite surface at different pressure
图11 不同压力下甲烷分子在孔径2 nm伊/蒙混层中沿z轴方向的密度剖面
Fig.11 Density profiles of methane molecule adsorption on 2 nm mixed-layer illite/smectite surface at different pressures
3 结 论
(1)构建了伊/蒙混层矿物的分子模型。利用分子模拟软件Materials studio模拟了甲烷在伊/蒙混层中的吸附行为，并与实验数据对比，证实了模型的可靠性。
(2)甲烷在不同孔径蒙脱石、伊/蒙混层和伊利石上的平均等量吸附热均小于42 kJ/mol，说明甲烷在黏土矿物表面发生物理吸附。甲烷在3种黏土矿物上的吸附等温线呈现出相同的规律，符合I型吸附曲线特征，Langmuir模型的拟合精度较高。3种黏土矿物对甲烷的吸附能力顺序为:蒙脱石>伊/蒙混层>伊利石。黏土矿物孔径越小，对甲烷的吸附能力越强，过剩吸附量越大，但受到空间限制，绝对吸附量及饱和吸附量随孔径减小而降低。
(3)在伊/蒙混层矿物中，层间离子比甲烷更靠近黏土片层，蒙脱石一侧吸附的甲烷密度大于伊利石一侧，表明蒙脱石片层的吸附能力强于伊利石。随着孔径增加，伊/蒙混层中的甲烷吸附量增加，吸附层数增加，吸附相密度减小。低压时，蒙脱石和伊利石片层上甲烷吸附密度的差异较为显著，但随着压力升高，吸附作用的影响减弱，导致两侧密度差异逐渐减小。
[1] TIAN H,LI T F,ZHANG T W,et al.Characterization of methane adsorption on overmature Lower Silurian-Upper Ordovician shales in Sichuan Basin,southwest China:Experimental results and geological implications[J].International Journal of Coal Geology,2016,156:36-49.
[2] CURTIS J B.Fractured shale-gas systems[J].AAPG Bulletin,2002,86(11):1921-1938.
[3] NUTTAL B C,CORTLAND E,BUSTIN M R,et al.Analysis of Devonian black shales in kentucky for potential carbon dioxide sequestration and enhanced natural gas production[J].Greenhouse Gas Control Technologies,2005,2(2):2225-2228.
[4] ROSS D J K,BUSTIN R M.The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs[J].Marine and Petroleum Geology,2009,26(6):916-927.
[5] LIU S G,Ma W X,JANSA L,et al.Characteristics of the shale gas reservoir rocks in the Lower Silurian Longmaxi Formation,East Sichuan Basin,China[J].Acta Petrologica Sinica,2011,31(2):187-220.
[6] ROSS D J K,BUSTIN R M.Characterizing the shale gas resource potential of Devonian-Mississippian strata in the Western Canada sedimentary basin:Application of an integrated formation evaluation[J].AAPG Bulletin,92,87-125.
[7] SCHETTLER J,PARMELY C,JUNIATA C.Contributions to total storage capacity in Devonian shales[A].SPE Eastern Regional Meeting[C].Lexington:1991:22-25.
[8] LU X C,LI F C,WATSON A T.Adsorption measurements in Devonian shales[J].Fuel,1995,74(4):599-603.
[9] ROSS D J K,BUSTIN R M.The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs[J].Marine & Petroleum Geology,2009,26(6):916-927.
[10] YANG F.Characterization of microscopic pore structures in shale reservoirs[J].Acta Petrolei Sinica,2013,34(2):301-311.
[11] JIAO S J，HAN H，WENG Q P,et al.Scanning electron microscope analysis of porosity in shale[J].Journal of Chinese Electron Microscopy Society,2012,31(5):432-436.
[12] KROOSS B M,BUSCH A,ALLES S,et al.Experimental investigation of molecular diffusion of CO2 in coals and shales[A].International Conference on Gas-Water-Rock Interactions Induced by Reservoir Exploitation,CO2 Sequestration,and other Geological Storage[C].Malmaison,2003.
[13] GODEC M,KOPERNA G,PETRUSAK R,et al.Potential for enhanced gas recovery and CO2 storage in the Marcellus Shale in the Eastern United States[J].International Journal of Coal Geology,2013,118:95-104.
[14] DAY-STIRRAT R J,LOCUKS R G,MILLIKEN K L,et al.Phyllosilicate orientation demonstrates early timing of compactional stabilization in calcite-cemented concretions in the Barnett Shale (Late Mississippian),Fort Worth Basin,Texas (U.S.A)[J].Sedimentary Geology,2008,208(1):27-35.
[15] YANG Feng,NING Zengfu,LIU H Q.Fractal characteristics of shales from a shale gas reservoir in the Sichuan Basin,China[J].Fuel,2014,115(1):378-384.
[16] 熊健，刘向君，梁利喜.甲烷在黏土矿物狭缝孔中吸附的分子模拟研究[J].煤炭学报,2017,42(4):959-968.
XIONG Jian,LIU Xiangjun,LIANG Lixi.Molecular simulation study on the adsorption behaviors of methane in slit-like clay mineral pore[J].Journal of China Coal Society,2017,42(4):959-968.
[17] 隋宏光，姚军.页岩黏土矿物CH4/CO2吸附规律的分子模拟[J].东北石油大学学报,2016,40(2):90-98,127.
SUI Hongguang,YAO Jun.Molecular simulation of CH4/CO2 adsorption in clay minerals[J].Journal of Northeast Petroleum University,2016,40(2):90-98,127.
[18] HAO Y Z,YUAN L F,LI P C,et al.Molecular simulations of methane adsorption behavior in illite nanopores considering basal and edge surfaces[J].Energy & Fuels,2018,32(4):4783-4796.
[19] JIN Z H,FIROOZABADI A.Methane and carbon dioxide adsorption in clay-like slit pores by Monte Carlo simulations[J].Fluid Phase Equilibria,2013,360:456-465.
[20] 李文华，房晓红，李彬，等.蒙脱石吸附CH4和CO2的分子模拟[J].东北石油大学学报,2014,38(3):25-30,6-7.
LI Wenhua,FANG Xiaohong,LI Bin,et al.Molecular simulation of the sorption of methane and carbon dioxide in the montmorillonite[J].Journal of Northeast Petroleum University,2014,38(3):25-30,6-7.
[21] 卢双舫，沈博健，许晨曦，等.利用GCMC分子模拟技术研究页岩气的吸附行为和机理[J].地球科学,2018,5(43):1783-1791.
LU Shuangfang,SHEN Bojian,XU Chenxi,et al.Study on adsorption behavior and mechanism of shale gas by using GCMC molecular simulation[J].Journal of Earth Science,2018,5(43):1783-1791.
[22] ZHAND J,CLENNELL M B,LIU K,et al.Methane and carbon dioxide adsorption on illite[J].Energy & Fuels,2016,30(12):10643-10652.
[23] SKIPPER N T,REFSON K,MCCONNELL J D C.Computer simulation of interlayer water in 2:1 clays[J].The Journal of Chemical Physics,1991,94(11):7434-7445.
[24] GEATCHES D L,WILCOX J.Ab initio investigations of dioctahedral interlayer-deficient mica:Modelling 1 M polymorphs of illite found within gas shale[J].European Journal of Mineralogy,2014,26(1):127-144.
[25] ALTANER S P.Comparison of structural models of mixed-layer illite/smectite and reaction mechanisms of smectite illitization[J].Clays & Clay Minerals,1997,45(4):517-533.
[26] CYGAN R T,LIANG J J,KALINICHEV A G.Molecular models of hydroxide,oxyhydroxide,and clay phases and the development of a general force field[J].The Journal of Physical Chemistry B,2004,108(4):1255-1266.
[27] MARTIN M G,SIEPMANN J L.Transferable potentials for phase equilibria.1.united-atom description of n-alkanes[J].The Journal of Physical Chemistry B,1998,102(14):2569-2577.
[28] INOUE A,VELDE B,MEUNIER A,et al.Mechanism of illite formation during smectite-to-illite conversion in a hydrothermal system Arsuvurr Iuoup[J].American Mineralogist,1988,73:1325-1334.
[29] JI L M,ZHANG T W,MILLIKEN K L,et al.Experimental investigation of main controls to methane adsorption in clay-rich rocks[J].Applied Geochemistry,2012,27(12):2533-2545.
[30] 安部郁夫,石川达雄，近藤精一.吸附科学[M].北京:化学工业出版,2006.
[31] ADAM N.Sorption and desorption of gases (CH4,CO2) on hard coal and active carbon at elevated pressures[J].Fuel,1998,77(11):1243-1246.
[32] EVERETT D H,POWL J C.Adsorption in slit-like cylindrical micropores in the Henry’s Law Region[J].Journal of the Chemical Society Faraday Transactions,1976,72(1):619-636.
[33] PENG D Y,ROBINSON D B.A new two-constant equation of state[J].Minerva Ginecologica,1976,12(11-12):3069-3078.
