页岩气多层吸附的分子模拟与理论模型

穆中奇，宁正福，吕方涛，顾可名，刘 蓓

(1.中国石油大学(北京) 油气资源与探测国家重点实验室，北京 102249； 2.中国石油大学(北京) 石油工程教育部重点实验室，北京 102249； 3.中国石油大学(北京) 重质油国家重点实验室，北京 102249)

引 言

目前，已经有许多模型用于描述页岩气的吸附现象，最常用的是Langmuir方程，该模型假设吸附质在均匀表面发生单层吸附[1]，但有学者认为该假设用在页岩吸附中不合理，因为页岩孔隙表面的性质是不均匀的[2]，且气体在微孔中发生填充型吸附[3]。即使用描述多层吸附的Brunauer-Emmet-Teller(BET)方程[4]，假设吸附表面均匀这一问题仍未解决。基于吸附势理论的Dubinin-Radushkevich(DR)和Dubinin-Astakhov(DA)方程[5]适用于描述各种方式的吸附及非均匀吸附，但它们不能给出吸附相的确切形态，难以用于吸附相流动问题的研究中。

利用分子模拟技术模拟页岩气在孔隙中的密度分布，可以直观地研究吸附方式[6-7]。许多页岩气分子模拟研究[8-10]表明，无论是矿物孔还是有机质孔，其中的吸附近似为单层吸附，且存在微弱的第二吸附层，特别是在微孔中，两侧孔壁作用叠加使得双层吸附更加明显。因此，为了更精细地表征页岩气吸附，需要探索如何将多层吸附模型应用到页岩气中。BET方程是经典的多层吸附模型，其最初是在亚临界条件下提出的，一些学者试图将其拓展到页岩气吸附这种超临界环境中。Wang等[11]将模型中的饱和蒸气压换成“饱和吸附压力”，但并未给出它的确定方法。Yu等[12]将饱和蒸气压换成“虚拟饱和蒸气压”，但它本身缺乏明确的物理意义。Zhou等[13]将模型中的相对压力p/p0换成游离气与吸附气的密度比ρg/ρa，但这基于吸附相密度恒定的假设，而这个假设却与分子模拟的结果[8-10]相悖，即吸附相密度应当是变化的。总之，目前需要更合理的超临界多层吸附模型。

为了明确页岩气的多层吸附特征并建立合适的理论模型，首先用分子模拟方法模拟了真实储层条件下有机质孔隙中甲烷的赋存状态，然后基于化学热力学理论提出了一种超临界气体多层吸附模型，并通过拟合分子模拟获得的吸附等温线验证模型的可靠性。

1 页岩气吸附的分子模拟

1.1 孔隙与吸附质模型

页岩气在有机质孔中赋存的分子模拟通过Materials Studio软件完成。以龙马溪组页岩为代表的四川盆地海相页岩的有机质主要为Ⅱ型并处于过成熟阶段。Ⅱ型干酪根本身氧/碳原子比很低，且随着其裂解和缩合反应不断变化，过成熟干酪根的氢/碳原子降低，芳构化程度升高，最终演化为石墨[14]。因此，石墨烯可以作为一种替代过成熟有机质的简化模型。熊健[15]提取了龙马溪组页岩干酪根并进行低温氮气吸附测试，吸附与脱附等温线形状表明干酪根孔隙表现出了狭缝孔特点。因此使用石墨烯组成的平行板孔代表有机质孔[9,16]。孔隙由两个平行的孔壁组成，每个孔壁包含四层石墨烯，相邻两层间距为0.34 nm，在x、y方向上的尺寸为4.92 nm×5.11 nm。孔宽定义为两孔壁距离最近的两层石墨烯的距离，有1 nm、2 nm、3 nm、5 nm四种。天然气以纯甲烷代表，甲烷使用全原子模型。孔隙与气体模型如图1所示。

图1 孔隙与气体模型Fig.1 Pore and gas model

1.2 模拟过程

根据四川盆地深层页岩气藏的地质条件[17]，将模拟的压力定为0～80 MPa，温度定为120 ℃和150 ℃，模拟3 nm孔中的页岩气吸附等温线，用于代表非微孔(孔宽大于2 nm)中的多层吸附情况。具体过程是用巨正则蒙特卡洛(GCMC)方法确定孔隙中甲烷分子数平均值，并用式

(1)

计算过剩吸附量。式中，nex为单位表面积过剩吸附量，mol/m2；S为吸附表面积，m2；NCH4为孔隙中甲烷分子数；MCH4为甲烷摩尔质量，kg/mol；NA为阿伏伽德罗常数；ρCH4为甲烷游离气密度，kg/m3；Vfree为孔隙自由空间体积，m3。自由空间体积定义为气体可充填的空间体积，可以用吸附实验中用到的氦气法测量[18]，这里用GCMC方法在30 ℃、1 MPa条件下测量，即

(2)

式中，NHe为孔隙中氦气分子数；MHe为氦气摩尔质量，kg/mol；ρHe为氦气密度，kg/m3。模拟与处理结果时甲烷与氦气的密度与逸度由美国国家标准与技术研究院(NIST)的REFPROP数据库计算。

另外，分别选择5 MPa、30 MPa、80 MPa代表页岩气开发过程中的低压、中压、高压状态。通过以下过程模拟120 ℃和150 ℃下上述几个压力状态的页岩气在1 nm、2 nm、3 nm、5 nm孔隙中的分布：首先用正则蒙特卡洛方法(CMC)将相应数目的甲烷分子填充入孔内(甲烷分子数根据GCMC的计算结果确定)，最终获得一个能反映平均吸附状态的能量最低构型，用于后续的动力学模拟；其次用分子动力学方法(MD)在正则(NVT)系综下优化孔内气体分布；最后通过统计不同时刻气体分子在孔内的位置获得密度分布剖面。

1.3 模拟细节

模拟时固定孔隙，允许甲烷自由运动，并使用周期性边界条件。模拟的时间步设为1 fs。蒙特卡洛模拟中，总步数设为5×106步，其中前2.5×106步用于体系平衡，后2.5×106步用于统计。MD模拟中，总步数设为1×106步，其中前2×105步用于体系平衡，后8×105步用于统计。用Nosé热浴法维持MD中体系温度稳定。力场选择COMPASS[19]。静电相互作用求和使用Ewald法，计算精度设为4.186 J/mol。范德华力求和使用Atom-based法，截断半径设为1.85 nm。

2 超临界气体多层吸附模型

亚临界条件下，吸附平衡时游离气的冷凝速率等于吸附气的蒸发速率，这可以作为推导吸附模型的基础。超临界条件下气体不会冷凝与蒸发，上述规律不能使用，但可以应用相平衡热力学理论。因此，本文从吸附相、游离相平衡时的热力学规律出发推导模型。

2.1 绝对吸附量模型

可以将多吸附层看作由气体分子在表面形成的大量多聚体组成[20]，如图2所示。定义Bj为由j个气体分子Bg组成的多聚体(j=1,2…,n)；n为最大吸附层数。

图2 吸附层的多聚体示意图Fig.2 Schematic diagram of multimers in adsorbed layer

真实气体的化学势为[21]

(3)

而大量Bj组成的吸附相成分的化学势为[20]

(4)

认为在多层吸附中存在以下平衡关系[20]：

BgB1；

(5)

Bg+BjBj+1。

(6)

式(5)表示形成第一层吸附层的平衡关系，式(6)表示形成第j+1层吸附层的平衡关系。

形成第一层吸附层且达到平衡时，摩尔吉布斯自由能变应为0，可以表达为

(ΔGm)1=μ1-μg=0。

(7)

将式(3)、(4)代入式(7)，并令j=1得到

(8)

定义标准摩尔吉布斯自由能变为

(9)

(10)

标准吸附平衡常数指当反应物与产物分别以各自标准态为参考态时的吸附平衡常数。对于一个确定的反应，标准平衡常数只与温度有关[21]。将式(9)、(10)代入式(8)可得到

(11)

形成第j+1层吸附层且达到平衡时，摩尔吉布斯自由能变应为0，且可以表达为

(ΔGm)j+1=μj+1-μj-μg=0。

(12)

将式(3)、(4)代入式(12)得到

(13)

定义标准摩尔吉布斯自由能变为

(14)

(15)

将式(14)、(15)代入式(13)得

(16)

由式(16)可推得θj、θj-1、……θ2的表达式，将他们代入式(16)可得到

(17)

(18)

(19)

θj+1=xjθ1(j=1,2,…)。

(20)

根据式(20)，各种多聚体的覆盖度之和为

(21)

总吸附量为

(22)

式中，nabs为单位表面积上的绝对吸附量，mol/m2；nm为单位表面积上单层最大吸附量，mol/m2。将式(21)代入式(20)得到

(23)

将式(23)代入式(22)得到超临界吸附模型，其与传统的BET模型[4]有相同的形式，即

(24)

(25)

计算。

2.2 过剩吸附量模型

分子模拟中，通过式(1)得到的吸附量均为过剩吸附量。孔内气体密度分布结果表明，页岩气在有机质孔中的吸附相体积保持恒定，因此可通过吸附相体积实现绝对吸附量与过剩吸附量之间的转换。可用式

(26)

(27)

计算。特别当n=1时，超临界多层吸附模型将退化为描述单层过剩吸附量的Langmuir模型[23]，即

(28)

3 结果与讨论

3.1 孔内气体分布

经MD优化后各温度下气体密度分布如图3所示。在非微孔中，从孔壁到孔中心，依次划分为以下几个区域(以宽为5 nm的孔为例在图3中标出)：①不可到达区。在不同孔宽、温度、压力下，不可到达区都没有甲烷分子，厚度均约0.18 nm，接近碳原子的范德华半径0.172 nm，这一区域是孔壁碳原子排斥产生的。②第一吸附层、第二吸附层。根据密度分布曲线出现的波动可以划分出两个吸附层。靠近孔壁的吸附层定义为第一吸附层，远离孔壁的吸附层定义为第二吸附层。通过在相应区域内求各个位置密度的算术平均值得到两吸附层的密度，与游离气密度理论值共同在图4中比较。可以看出，第一吸附层的密度远高于第二吸附层，第二吸附层的密度只比游离相稍高一些，这说明第一吸附层受到的吸附作用很强，而第二吸附层受到的吸附作用很弱。在不同孔宽、温度、压力下，两吸附层均各厚约0.38 nm，等于甲烷分子的直径。3 nm孔中两吸附层的密度基本等于5 nm孔而略低于2 nm孔，说明孔宽超过2～3 nm，两侧孔壁对气体的叠加作用力就很弱了，两吸附层密度将不受孔宽影响，因此，3 nm孔的吸附状况可以代表非微孔中的吸附状况。③游离相。这一区域的密度剖面稳定在游离气密度理论值上。孔宽越大，游离相范围也越大。在孔宽为1 nm的孔中，受空间所限，只存在不可到达区和填充了全部自由空间的吸附相，分别厚约0.18 nm和0.64 nm，且厚度不随温度、压力变化。因此，页岩气在有机质孔中的吸附相体积保持不变，在非微孔中发生双层吸附。

图3 甲烷在有机质孔中的密度剖面与游离相密度理论值Fig.3 Density profiles of methane in organic pores and theoretical density of free gas

图4 非微孔中吸附层密度与游离相密度Fig.4 Density of adsorbed layers and free gas in non-micropores

密度剖面的形状受压力影响较大。在非微孔中，压力从5 MPa升到30 MPa再升到80 MPa，第一、第二吸附层的分界点从一个骤然增长的拐点变成了波谷，又变成了更加明显的波谷；而第二吸附层与游离相之间的分界点从一个缓慢增长的拐点变成了一个骤然增长的拐点，又变成了一个轻微的波谷；第一吸附层内始终有一个明显的波峰；第二吸附层内从没有波峰变成了出现波峰，又变成了更加明显的波峰。在宽1 nm的孔中，吸附层内的波峰、吸附层之间的波谷随压力升高越来越明显。总体来说，随着压力升高，吸附层内波峰不断强化，不同区域分界点处波谷逐渐出现，这表明吸附层内的分子排列得愈加整齐有序，吸附层间及其与游离相间的界限愈加分明，孔壁附近气体分布的连续性变差，这可以归因于吸附层密度升高会加强对自身及周边分子的控制与吸引。

3.2 多层吸附等温线

页岩气在3 nm有机质孔中的多层吸附等温线如图5所示。随着压力增大，过剩吸附量早期增长，中后期下降，这是不同压力阶段吸附相密度与游离相密度增大规律不同造成的。在过剩吸附量下降段，其下降趋势会随压力升高而减缓，这一点与真实页岩的实测高压(超过50 MPa)吸附等温线一致[24-25]。随着温度升高，过剩吸附量降低，这表明页岩气在有机质孔中的吸附是放热过程，温度升高能促进吸附气解吸。

图5 分子模拟过剩吸附等温线与模型拟合曲线Fig.5 Simulated excess adsorption isotherms and model fitting isotherms

3.3 模型拟合结果

表1 吸附等温线的参数拟合结果Tab.1 Fitting parameters of adsorption isotherms

用超临界多层吸附模型计算的两吸附层密度随压力的变化趋势与分子模拟结果大致吻合，如图6所示。但可以看出密度预测值存在一定误差，表现为：第一吸附层的模型计算值偏大，第二吸附层的模型计算值偏小；低压下误差较大，高压下误差较小。实际上，无论是Langmuir模型还是多层吸附模型，它们所基于吸附相状态方程都简化了吸附层内分子间及吸附层间的相互作用[18]，而在分子模拟及真实吸附中这些作用将影响吸附层密度，这导致了模型预测的吸附层特征不够准确。

图6 分子模拟与超临界多层吸附模型计算的吸附层密度比较Fig.7 Comparison of adsorption layer density obtained by simulation and by supercritical multilayer adsorption model

4 结 论

(1)页岩气在有机质非微孔中为双层吸附，第一吸附层受到的吸附作用远强于第二吸附层；每个吸附层的体积不随压力、温度和孔宽变化；当孔宽超过2～3 nm，吸附层密度不受孔宽影响。在宽1 nm的有机质孔中，仅存在填充整个孔隙空间的吸附相，不随压力和温度变化。随着压力升高，吸附层内的分子排列得更加整齐，吸附层间及其与游离相间的界限更加分明。

(2)深部储层条件下，随着压力增大，页岩气在有机质非微孔中的过剩吸附量先增长、后下降，且下降趋势会随压力升高而减缓。

(3)当最大吸附层数为两层时，超临界多层吸附模型的拟合结果与Langmuir方程同样良好，但得到的吸附层厚度比Langmuir方程更准确，且计算出的各吸附层密度随压力变化的趋势更准确。