引言

近年来，页岩气在世界范围内受到广泛关注，美国和加拿大已经实现了页岩气的商业开采^[1]并取得较好成果，中国在页岩气的研究现状、成藏机理和有利区的评价方面做了大量工作，但近年才开始研究页岩气的吸附现象^[2-3]。页岩气赋存状态包括游离态、吸附态和溶解态^[4]，主要由甲烷气体组成，含量甚至可达95%以上^[5]，大量学者利用甲烷来近似代替页岩气并得到广泛应用^[6]。页岩岩石物理特性一直是研究热点^[7-8]，其矿物组成以石英为主^[9]。不同储层的吸附性能受环境因素以及页岩本身物理化学性质影响，其中，润湿性是岩石物理特性的一个重要参数。润湿性是指不相混的两相流体与岩石固相表面接触时，其中某一相流体沿着岩石表面铺开的现象^[10]。由于储层岩石的润湿性非常复杂，一般采用室内实验测试来确定岩石的润湿性。Bai等运用接触角法测量了Fayetteville页岩气储层润湿性，结果表明该页岩储层最初表现为中性润湿^[11]。Lakatos等运用接触角法测量了非常规气藏储层的润湿性^[12]。刘向君等使用接触角法研究分析了龙马溪组页岩润湿性，研究结果表明龙马溪组页岩表面既亲油又亲水，孔隙表面处润湿性存在差异，出现斑状润湿^[13]。

随着研究的深入，学者们除了从宏观实验方面^[14-16]研究气体与页岩矿物之间的相互作用外，还进行了微观层面的研究，而分子模拟是一种将理论方法与计算技术相结合来模拟分子微观运动行为的研究方法，并已成为研究流体在固体表面的有效手段。研究者们采用分子模拟建立矿物结构并进行了不同材料表面润湿性的模拟研究^[17-19]。Chai等采用分子动力学模拟研究了水在不同羟基和烷基化的无定形二氧化硅表面上的微观润湿性^[20]，研究结果表明，疏水性随着表面硅烷化程度增大而增大，表面硅烷化的比例超过70%的无定形二氧化硅表面接触角基本保持不变。Iglauer等使用分子动力学研究了盐水-二氧化碳的界面张力和盐水-二氧化碳-岩石表面水接触角^[21]。刘冰等采用分子动力学方法研究了NaCl对水滴在砂岩表面的微观润湿行为的影响^[22]。以上采用分子模拟方法进行的研究，为表面润湿性修饰条件下的模拟提供了借鉴。基于前人的研究，选择甲烷气体模拟页岩气，并以石英为基本矿物载体，通过润湿性修饰石英表面构建不同润湿性模型，采用分子模拟方法研究甲烷在不同润湿性石英表面的吸附行为，并通过计算等温吸附线、密度分布分析甲烷的微观吸附机制。

1 模型与方法 1.1 模型的建立

采用Material Studio 2017软件数据库进行甲烷与矿物模型建立，石英原胞晶格参数为：$a$=0.491 3 nm，$b$=0.491 3 nm，$c$=0.540 5 nm，$\alpha$=90°，$\beta$=90°，$\gamma$=120°，结构如图 1所示。

文献中^[23]通常直接对石英原胞超胞化后建层，未考虑不同晶面的影响。在构建石英超胞时，观察发现截取不同晶面存在一定差异，以2$\times$2$\times$2的超胞为例，截取(100)(010)(001)晶面均需断开8个Si—O键，如图 2所示，但不同晶面提供的吸附位具有差异。

分别对原始晶胞截取两层厚度的(100)，(010)和(001)面，建立5$a{\times}5b$的超晶胞结构，构建无官能团修饰狭缝模型，其中，狭缝宽度分别设置为1 nm与2 nm，如图 3所示。

对不同晶面构建模型进行甲烷吸附模拟，对比得到吸附量最大的模型后确定该晶面，再对原始晶胞截取该晶面，建立6$a{\times}6b$的超晶胞结构，并对晶面进行羟基化和甲基化修饰^[24-25]，构建亲油(甲基)和亲水(羟基)模型。

1.2 研究方法

蒙特卡罗方法(GCMC)广泛应用于吸附研究，文中的吸附计算在Materials Studio 2017的Sorption模块中进行，采样法为Metropolis，使用Fixed pressure(恒温定压逐点计算)任务来进行模拟，选取Universal力场，库伦力相互作用采用Ewald法，范德华力相互作用采用Atom Based法，且L--J势能截断距离为1.25 nm，模拟的总步数为3$\times$10$^{{\rm 6}}$步，平衡步数为1.5$\times$10$^{{\rm 6}}$步，前一半模拟步数用于保证吸附模拟达到吸附平衡，后一半模拟步数用于进行系综平均，具体参数设置与文献[23]相同。

借鉴实际页岩储层温度、压力^[26]，拟定模拟温度为333 K，同时设置313 K进行不同温度对比，压力范围为0.1$\sim$40.0 MPa。

气体分子的逸度根据Peng-Robin状态方程^[27]进行计算。

吸附模拟后根据分子动力学基本原理，采用Forcite模块下的Dynamics进行动力学模拟，模拟过程在三维周期性条件下进行，以减小模拟体系尺寸效应。

分子的初始速度由Maxwell-Boltzmann分布随机产生，选用NVT系综，时间步长为1 fs，总步长为2 ns，前1.0 ns用来平衡模拟体系，后1.0 ns用来采集数据进行分析^[23]，控温方法为Nose，模拟结束后进行时平均密度分析。

2 结果与讨论 2.1 不同晶面模型吸附量

不同温度和压力下，甲烷在不同晶面狭缝模型中的吸附等温线如图 4所示。

由图 4可知，甲烷的吸附量随着压力的升高而增大，在低压下的吸附速度更快，高压下的吸附速度逐渐减慢。1 nm模型的吸附量趋于饱和，而2 nm模型的吸附量呈继续上升的趋势。吸附量随温度的升高而减少，且温度对吸附量的影响程度随着孔径增大而增大；大孔径模型的吸附量远大于小孔径模型的吸附量。在温度、压力固定的条件下，(100)面构建的2 nm模型中甲烷的吸附量大于相同孔径的其他两个模型的吸附量，且不同温度下吸附量大小依次是(100)>(010)>(001)模型；构建的1 nm模型中甲烷的吸附量在313 K条件下依次也是(100)>(010)>(001)模型，而333 K条件下是(100)≈(010)>(001)，因此，研究使用石英模拟页岩主要矿物进行吸附研究时，不能单纯地直接对原胞超胞化后就建层，还应该进行不同截面分析，考虑不同面优化后的稳定性与能量大小，再进行后续吸附等研究。

2.2 润湿性修饰模型中甲烷吸附量

基于不同晶面模型吸附量结果，对(100)面底端氧原子层进行甲基和羟基等不同润湿性修饰，得到图 5所示不同润湿性修饰模型，1 nm模型盒子大小为2.456 5 nm$\times$2.702 6 nm$\times$3.051 7 nm (羟基化)，2.456 5 nm$\times$2.702 6 nm$\times$3.389 5 nm(甲基化)；2 nm模型盒子大小为2.456 5 nm$\times$2.702 6 nm$\times$4.051 7 nm (羟基化)，2.456 5 nm$\times$2.702 6 nm$\times$4.389 5 nm(甲基化)。

图 6为不同润湿性修饰的石英模型中甲烷的吸附等温线。

从图 6可以看出，在1 nm孔模型中，无论温度为313 K还是333 K，0.1 MPa到10.0 MPa压力段甲烷的吸附量均快速增加，10.0 MPa后增长速度变缓并趋于饱和，甲基修饰表面构建的模型中吸附量大于羟基修饰表面构建的模型中甲烷吸附量，表明甲基修饰表面更加亲甲烷。在2 nm模型中，0.1 MPa到20.0 MPa压力段甲烷吸附量快速增加，20.0 MPa后增长速度变缓，不同于1 nm孔隙模型的是增长趋势只是趋缓而未出现饱和，甲基修饰表面构建的模型中吸附量也略大于羟基修饰表面构建的模型中甲烷吸附量。在一定压力条件下，不论是甲基还是羟基修饰的表面模型中，吸附量均随温度的升高而减小，且大孔径模型中温度的影响更加显著。此外，与相同条件下无润湿性修饰模型中的吸附量对比可知，润湿性修饰后的吸附量均远大于未修饰模型中的吸附量，羟基修饰模型中的吸附量增大源于官能团增强了与甲烷之间的相互作用，且模拟中未加入水分子进行竞争吸附。

2.3 润湿性修饰模型中甲烷密度分布

313 K条件下不同模型中甲烷密度分布如图 7所示。由图 7可知，甲烷在修饰的1 nm模型中的密度分布存在略微差异，压力为5.0 MPa时，羟基修饰模型中的甲烷在壁面处形成单吸附层，随着压力增加，表层吸附相甲烷逐渐趋于饱和。甲基修饰模型中密度分布规律类似，在5.0 MPa后，甲烷同样在壁面处形成单吸附层，且密度大于相同条件下羟基模型中的密度，表明甲基修饰模型吸附能力更强，当压力大于10.0 MPa后，两种模型中的密度均有缓慢增加。不同官能团修饰的2 nm模型中甲烷密度分布趋势一致，压力为5.0 MPa时，甲烷在壁面处形成单吸附层且甲基修饰模型中的密度更大；5.0 MPa后随着压力增加，密度分布峰值逐渐增大，吸附层数逐渐增多，孔隙中甲烷密度分布有两个主要峰值和3个次峰，近壁面处出现明显的吸附层，孔隙中央的吸附层明显弱于壁面处吸附层；高压下壁面处的吸附层趋于饱和，孔隙中的吸附层更趋明显，密度分布峰值显示各模型中均形成了5个吸附层。

3 结论

（1）不同晶面模型中的吸附等温线表明(100)面模型的吸附量大于其他两个模型的吸附量，使用石英进行研究时，不能直接对原胞进行超胞化后就建层，还应该进行不同截面分析研究。

（2）甲烷在不同润湿性修饰模型中的吸附量均随着压力的增加而增加，吸附量随温度的升高而减小，温度对甲烷的吸附量的影响小于孔径对甲烷的影响。相同条件下，润湿性修饰后模型的吸附量均远大于未修饰模型中的吸附量。

（3）润湿性修饰1 nm模型中的密度分布表明甲烷只在壁面处形成单吸附层，而润湿性修饰2 nm模型中的甲烷形成多个吸附层。