深层页岩伊利石中甲烷吸附特征分子模拟

黄亮，陈秋桔，吴建发，杨琴，张鉴，黄山，周文，邹杰

(1.油气藏地质及开发工程国家重点实验室(成都理工大学)，四川 成都，610059；2.成都理工大学 能源学院，四川 成都，610059；3.中石油西南油气田分公司 页岩气研究院，四川 成都，610051)

经过近10年的努力探索，川南地区3 500 m以浅页岩气已实现规模效益开发，页岩气开发逐渐向埋深处于3 500～4 500 m 的深层页岩气迈进，该埋深范围的页岩气资源量高达8.7×1012m3，占埋深4 500 m 以浅资源总量的87%，是下一步勘探开发的重点领域[1]。目前深层页岩气开发理论尚不完善，开发技术仍欠成熟，为实现其规模效益开发，迫切需要强化相关基础理论的科研攻关[2]。

页岩气在储层中以吸附、游离和溶解3种状态赋存，其中吸附气含量可占气体总量的20%～85%[3]。黏土矿物是页岩气主要的赋存空间之一，其中，伊利石对甲烷吸附具有重要贡献。唐鑫等[4]发现四川盆地龙马溪组页岩中黏土矿物质量分数可达到42%，而黏土矿物中伊利石质量分数可高达84%；齐荣荣[5]通过低温二氧化碳吸附表征伊利石矿物微孔结构，发现伊利石中存在丰富的微孔，微孔体积可超过总孔体积的10%。因此，明确页岩伊利石中甲烷吸附特征对页岩气藏储量评估具有重要意义。

室内吸附实验是研究页岩气吸附特征的直接手段之一，大量学者通过吸附实验测试甲烷在页岩中的吸附特征[6-10]。深层页岩气储层温度可超过150 ℃，压力可超过90 MPa。虽然目前部分吸附实验装置的测试温度可达到150 ℃，但实验装置的最大测试压力尚无法达到实际深层页岩气储层条件。

分子模拟方法可以突破吸附实验装置的温压耐受极限，能够从分子角度认识和解释气体吸附机理，研究页岩纳米级孔隙中气体吸附行为。针对甲烷在页岩岩石组分中的吸附特征，国内外学者采用分子模拟方法开展了大量研究[11-14]。但存在以下不足：

1) 黄亮等[11,13-14]主要针对3 500 m 以浅页岩气储层，模拟的温度和压力均较低；

2) 杨建等[12]在深层高温、高压条件下开展模拟研究，但由于采用巨正则蒙特卡洛(GCMC)算法，未考虑受孔喉连通程度影响的分子热运动行为，与分子实际运动规律不符；

3) 大部分模拟研究[11-14]基于气体状态方程，结合压力和逸度系数，近似获取气体的化学势，与实际吸附模拟中，由经验力场确定的气体化学势存在不一致性。这种不一致性导致的化学势偏差随压力增加而增大[15]，无法应用于深层高温高压条件。

4) HO等[16]通过Widom方法模拟化学势和压力的关系，主要适用于低密度条件，但不适用于深层页岩气储层。

准确计算高温高压条件下气体化学势难度很大。为通过分子模拟方法研究甲烷在深层页岩伊利石中的吸附特征，需首先获取准确的甲烷化学势。为此，本文首先采用GCMC 和分子动力学(MD)耦合的方法(GCMC-MD 方法)，开展不同温度和压力下甲烷在体相系统中的吸附模拟研究，构建一套甲烷化学势图版；其次，结合甲烷化学势与温度和压力的关系，建立一套甲烷化学势预测模型，结合确定的甲烷化学势，模拟甲烷在伊利石纳米孔隙中的吸附行为，分析不同温度和孔径下的甲烷吸附曲线特征；最后，基于模拟结果，分析压力、温度和孔径对伊利石中甲烷微观吸附特征的影响规律。

1 分子模拟方法

1.1 伊利石纳米孔隙分子模型

以川南泸州区块深层龙马溪组页岩气储层为研究对象，目标页岩矿物分析表明，伊利石是主要的无机黏土矿物。本文选择页岩气储层中伊利石转化过程末期的伊利石类型开展模拟，其单胞结构由GEATCHES等[17]构建。图1所示为构建的伊利石纳米孔隙分子模型，其单胞参数为a=5.22×10-10m，b=9.00×10-10m，c=10.32×10-10m，α=90.26°，β=103.05°，γ=89.97°，分子式为K(Si7Al)Al4O20(OH4)。伊利石中存在类质同象置换现象，一侧四面体上每8个Si原子中有1个被Al原子取代，导致伊利石结构呈电负性，电负性由层间的K+中和。

为避免破坏伊利石晶体结构，沿伊利石层状结构建立狭缝孔隙。首先，基于伊利石单胞建立超晶胞，超晶胞由8个×4个×2个单胞构成；然后，沿超晶胞(001)解理面切面形成狭缝孔壁；最后，将2 层孔壁组合形成伊利石狭缝孔隙，2 层孔壁表面氧原子中心平面之间的距离定义为狭缝孔径。图1所示为构建的伊利石狭缝孔隙分子模型，K+吸附于狭缝表面，可降低狭缝的有效孔隙尺寸。

1.2 甲烷吸附模拟细节

基于LAMMPS 分子模拟程序，采用GCMCMD 方法模拟甲烷在伊利石中的吸附。其中，GCMC 算法实现甲烷分子在伊利石中的插入和删除，MD算法实现甲烷分子在伊利石中的分子热运动。伊利石模型由ClayFF 力场描述，甲烷采用TraPPE 联合原子模型。在吸附模拟中，伊利石孔壁结构保持固定，狭缝表面K+可自由运动。

吸附模拟中以化学势作为输入参数，从而控制体系压力。化学势与压力的关系利用GCMCMD 方法在体相体系中获取。基于GCMC-MD 方法，结合确定的甲烷化学势，模拟伊利石孔隙中甲烷吸附行为。GCMC-MD 耦合模拟包括4 000 个循环，每个循环由初期的10 000个MD步和后期的2 500 个GCMC 步组成。其中，MD 步在正则系综下模拟，温度由Nosé-Hoover方法控制，截断半径为14×10-10m，时间步长为1×10-12s，体系采用周期性边界条件。吸附模拟得到伊利石中甲烷总气体量，其在平衡阶段的标准相对误差小于5%。过剩吸附量由总气体量换算得到。

式中：ne为过剩吸附量；nt为总气体量；ρ为甲烷体相密度；Vp为孔隙中甲烷可接触体积，由Metropolis蒙特卡洛积分法计算。

2 甲烷吸附化学势确定

2.1 甲烷化学势图版

图2所示为获取的甲烷化学势图版，温度范围为30～150 ℃，压力范围为0～110 MPa。图2中，a1，b1，c1和d1为待确定的回归常数；μ为化学势，kJ/mol；P为压力，MPa；T为温度，℃。由图2(a)可见：在不同温度下，甲烷体系压力与化学势均呈指数函数关系，随化学势增加，压力呈先缓慢增加后急剧增加的趋势。由图2(b)可见：在不同压力下，甲烷体系温度与化学势均呈线性函数关系，随化学势增加，甲烷体系温度逐渐降低。

图2 甲烷化学势图版Fig.2 Type curves of methane chemical potential

图3所示为甲烷化学势图版的验证结果。由图3(a)可见：甲烷体系密度随压力增加而逐渐增大，随温度升高而逐渐降低。在不同温度下，甲烷体系密度模拟结果与美国国家标准与技术研究院(NIST)数据库结果[18]基本一致。该结果表明：模拟的甲烷体系中，化学势-压力-密度具有准确的对应关系，构建的甲烷化学势图版准确和合理。

图3 甲烷化学势图版验证结果Fig.3 Validation results of type curves of methane chemical potential

图2中构建的甲烷化学势图版是基于边长为45×10-10m的立方体盒子模拟得到。为了检验化学势图版是否具有尺寸效应，基于边长为90×10-10m的立方形盒子重复甲烷吸附模拟研究，模拟温度为150 ℃，压力超过100 MPa。由图3(b)可见：随着模拟进行，甲烷体系的平均压力和平均密度逐渐收敛。不同模型尺寸下的体系平均压力和平均密度基本一致，表明构建的甲烷化学势图版无尺寸效应，具有代表性，适用于不同尺寸的甲烷体系。

2.2 甲烷化学势预测模型

图2表明甲烷体系压力与化学势呈指数函数关系；甲烷化学势与体系压力呈对数关系；在不同压力条件下，甲烷化学势与温度呈线性关系。基于化学势分别与压力和温度的关系，建立化学势与压力对数及温度的多元线性回归关系式：

式中：e1，f1和g1为待确定的回归常数。

为构建化学势与压力和温度的多元回归模型，将模拟的化学势数据分成回归组和验证组，其中回归组包括60，90，120，150 ℃这4 个温度下的化学势，验证组为30 ℃下的化学势。将低温条件下的化学势作为验证组，可检验回归模型的外延预测能力。基于回归组化学势，结合式(2)开展多元回归，建立甲烷化学势与压力和温度的定量关系模型：

基于式(3)预测不同温度和压力条件下的甲烷化学势，并与分子模拟结果对比，结果如图4所示。由图4可见：针对回归组数据，甲烷化学势的预测结果和模拟结果基本重合；针对验证组数据，甲烷化学势的预测结果和模拟结果吻合较好，相对误差低于7%，表明建立的甲烷化学势预测模型具有较高的预测精度。

图4 甲烷化学势模拟结果(实线)与预测结果(虚线)对比Fig.4 Comparison of methane chemical potential between simulation results (solid lines) and prediction results (dashed lines)

3 甲烷吸附曲线特征

3.1 压力和温度影响规律

基于确定的甲烷化学势，开展甲烷在伊利石纳米孔隙中的吸附模拟。为反映深层页岩气的吸附特征，模拟温度最高为150 ℃，压力最高超过90 MPa。图5所示为不同温度条件下，甲烷在孔径为2 nm 时的伊利石孔隙中的吸附曲线。甲烷总气体量包括孔隙表面吸附气量以及孔隙中间游离气量。随压力增加，甲烷总气体量呈初期迅速增加而后期缓慢增加的趋势，深层页岩储层高压条件下，纳米孔隙中甲烷总气体量未达饱和。这是因为在低压条件下，吸附气量和游离气量均随压力而增加；而在高压下，虽然吸附气量逐渐达到饱和，但游离气量随压力继续增加。过剩吸附量随压力增加呈先增大后减小的趋势，过剩吸附量在20 MPa 附近达到最大值，在80 MPa 后降为负值。根据过剩吸附量定义，过剩吸附量下降是由于甲烷体相密度随压力的增加速度大于吸附相密度的增加速度；而当体相密度大于吸附相密度时，过剩吸附量变为负值[19]。

图5 不同温度下甲烷吸附曲线(孔径为2 nm)Fig.5 Methane adsorption isotherms under various temperatures(pores with diameter of 2 nm)

在研究压力范围内，甲烷总气体量随温度增加而单调降低。在压力低于40 MPa 时，甲烷过剩吸附量随温度增加而单调降低；当压力高于40 MPa 时，不同温度下的甲烷过剩吸附量无明显差异。随着温度减小，最大过剩吸附量对应的压力略有减小。因此，相比中深层(埋深小于3 500 m)页岩气储层(较低温度和压力条件)，深层(埋深大于3 500 m)页岩气储层对甲烷的吸附能力大幅度降低。唐书恒等[20]发现伊利石中甲烷的过剩吸附量在20 MPa 下未达饱和，且随温度升高而降低。本文在对应温度压力范围的模拟结果与唐书恒等[20]的实验结果一致，进一步验证了模拟的合理性。

3.2 孔径影响规律

页岩气储层纳米孔隙十分发育，是气体主要的储集空间[21-22]。部分页岩的平均孔径仅为3.5 nm[15]。图6所示为不同伊利石孔径下的甲烷吸附曲线(120 ℃)。为消除比面积的影响，不同孔径模型中比表面积均相同。在研究压力范围内，甲烷总气体量随孔径增加而单调增加。在孔径分别为2 nm和4 nm孔隙内，甲烷过剩吸附曲线基本重合。当压力小于10 MPa 时，孔径分别为1，2 和4 nm 的孔隙内的过剩吸附量接近；当压力大于10 MPa时，孔径为1 nm的孔隙内的过剩吸附量低于孔径为2 nm 和4 nm 孔隙内的结果。在孔径为1 nm孔隙中，甲烷过剩吸附量约在10 MPa下达到最大；在孔径为2 nm和4 nm孔隙内，甲烷过剩吸附量约在20 MPa 下达到最大。伊利石孔径对甲烷过剩吸附的影响在小孔径范围更为明显。

图6 不同伊利石孔径下甲烷吸附曲线(120 ℃)Fig.6 Methane adsorption isotherms in illite with various pore sizes(120 ℃)

4 甲烷微观吸附特征

以孔径为4 nm 伊利石孔隙模型为例，研究甲烷在高温(120 ℃)高压(89.2 MPa)条件下吸附平衡后的微观分布特征，如图7所示。甲烷分子在狭缝孔隙的2个壁面附近各形成1个吸附层，吸附层内甲烷密度远大于孔隙中间密度。在狭缝孔隙中间，甲烷分子以游离的体相状态存在，不受孔隙壁面的影响。模拟的甲烷体相密度与NIST 数据库结果[18]基本一致，表明模拟结果合理和准确。

此外，吸附相和游离相之间存在1 个过渡区，在密度剖面上呈1个小凹凸结构，凹陷区域密度低于体相密度，凸起区域密度略高于体相密度，但远低于吸附层密度。过渡区小凹凸结构的形成包括3个原因。

1) 孔壁对过渡区甲烷分子的吸引力造成了小凸起结构。由于吸附层的空间位阻效应，小凸起区域甲烷分子与孔隙壁面的相互作用减弱，使得甲烷的密度大幅降低。

2) 吸附层分子与过渡区甲烷分子之间存在的排斥力造成了小凹陷结构。分子间存在非键相互作用，随分子间距离降低，相互吸引变为相互排斥，这也是壁面上吸附层密度呈逐渐增加趋势的原因之一。

3) 吸附层分子的不规则排列加剧了小凹陷结构的形成。由图7(b)可见：由于分子热运动，吸附层内分子并不呈直线规则排列。弯曲的吸附层使得其邻近区域的甲烷分子在凹凸的粗糙表面聚集，从而导致凹陷区单位厚度内的甲烷密度低于体相密度。

图7 伊利石纳米孔隙中甲烷微观分布(4 nm，120 ℃，89.2 MPa)Fig.7 Microscopic distribution of methane in illite nanopore(4 nm,120 ℃,89.2 MPa)

4.1 压力和温度影响规律

图8所示为不同压力和温度下的甲烷密度剖面和吸附层峰密度。由图8(a)可见：在1.8 MPa 和10.0 MPa 下，甲烷分子在2 个狭缝壁面上各形成1个吸附层。随压力继续增加，甲烷分子在壁面形成第2个吸附层。吸附层内的甲烷密度随压力增加而逐渐增大。由图8(c)可见：狭缝壁面上的第2个吸附层在压力大于10.0 MPa后形成。2个吸附层内的吸附层峰密度随压力增加呈相似的增加趋势，吸附层峰密度在高压下(89.2 MPa)增加幅度变小，呈现逐渐饱和的趋势。第1个吸附层内的吸附层峰密度始终大于第2个吸附层内结果。

图8 不同压力和温度下甲烷密度剖面和吸附层峰密度(2 nm)Fig.8 Methane density profiles and peak densities of adsorption phase under various pressures and temperatures(2 nm)

由图8(b)可见：在不同温度的高压条件下(约90 MPa 时)，甲烷分子在孔径为2 nm 孔隙的每个壁面附近形成2个吸附层。由于孔径较小，孔隙内未出现明显的体相区域。随温度增加，孔隙内各位置处的甲烷密度逐渐减小。由图8(d)可见：甲烷吸附层峰密度随温度增加而减小，其在低温范围内降低幅度较小，在高温条件下呈线性下降趋势。

4.2 孔径影响规律

图9所示为不同孔隙孔径下的甲烷密度剖面和吸附层参数。由图9(a)可见：在孔径为1 nm 的孔隙中，2个狭缝壁面对孔隙内甲烷分子存在强烈的耦合作用，孔隙内仅存在1个吸附层，吸附层的厚度约为0.5 nm。这是由于孔隙的孔径较小，狭缝壁面之间无足够空间形成2个吸附层，造成2个壁面的吸附层融合为1 个吸附层。甲烷在孔径为1 nm的孔隙中的吸附机理为微孔填充。在孔径为2 nm孔隙中，每个孔隙壁面附近存在2个吸附层，靠近壁面的吸附层峰密度略大于孔径为1 nm 的孔隙中的吸附层峰密度。由于2个狭缝壁面对孔隙内甲烷分子仍存在耦合作用，每个壁面附近形成第2个吸附层，其密度高于甲烷体相密度。

由图9(b)可见：在介孔中，随着孔径增加，2个壁面对孔隙内甲烷分子的耦合作用变弱，孔径为2 nm 的孔隙壁面上的第2 个吸附层(图9(a))退化为介孔中的过渡区小凹凸结构(图7(a))，介孔中每个壁面上均为单层吸附。随介孔的孔径增加，孔隙中间的甲烷游离相密度基本不变，与体相条件下的甲烷密度相同，介孔孔径增加主要增大孔隙内游离相的体积。

图9 不同孔径下甲烷密度剖面和吸附层参数(120 ℃，89.2 MPa)Fig.9 Methane density profiles and adsorbed layer parameters under various pore sizes(120 ℃,89.2 MPa)

由图9(c)可见：在孔径较小的孔隙中，狭缝壁面的耦合作用可影响甲烷的吸附层峰密度和吸附层厚度，吸附层峰密度与吸附层厚度之间存在反相关性。在孔径较大的介孔中，吸附层峰密度和吸附层厚度均不随孔径变化，吸附层厚度与甲烷分子动力学直径(0.38 nm)接近，呈典型的单层吸附特征。

5 结论

1) 建立了一套甲烷化学势与压力和温度的关系图版，甲烷化学势图版无尺寸效应，化学势-压力-密度具有准确的对应关系；构建了一套甲烷化学势预测模型，预测相对误差低于7%，具有较强的外延预测能力。

2) 在深层页岩气储层条件下，纳米孔隙内总气体量未达饱和，低压下以吸附气为主，高压下以游离气为主；相比中深层页岩气储层，深层页岩气储层过剩吸附气含量更低，游离气含量更高；在高压条件下(大于40 MPa)，温度对甲烷过剩吸附无明显影响。

3) 孔径对甲烷过剩吸附的影响在小孔径范围及高压条件下更为明显。在低压条件下，不同孔径中的甲烷过剩吸附量较为接近；当孔径大于4 nm 时，介孔孔径对甲烷吸附无明显影响。在深层页岩气储层条件下，甲烷在孔径为1 nm 的孔隙中仅存在1个吸附层，吸附机理为微孔填充；在孔径为2 nm 的孔隙的每个壁面上存在2 个吸附层；在孔径为4 nm 及更大孔径的介孔中，甲烷在每个壁面上单层吸附。

4) 在深层页岩气储层条件下，介孔中吸附相和游离相之间存在过渡区，过渡区小凹凸结构的形成包括3个原因，即孔壁对过渡区甲烷分子的吸引力、吸附层分子与过渡区甲烷分子之间的排斥力以及吸附层分子的不规则排列。