页岩有机质——石英复合孔隙中甲烷吸附与流动行为的分子模拟研究

Molecular Simulation of Methane Adsorption and Transport in Shale Organic Matter—Quartz Composite Nanopores

摘要: 吸附气是页岩气赋存在页岩储层中的主要方式之一，而且吸附气是页岩气后期产量的主要来源，页岩储层甲烷高温高压条件下的赋存特征是准确评估页岩气储量的关键。依据四川盆地五峰组——龙一₁亚段页岩无机地球化学特征，构建了石英、干酪根和复合孔组成的有机质——石英复合孔隙体系分子模型。通过巨正则蒙特卡洛(GCMC)、分子动力学(MD)和非平衡分子动力学(EFMD)模拟方法，研究了甲烷分子在复合孔隙体系中的微观赋存特征及其流动扩散特征。研究结果显示，随着压力的升高，赋存于孔隙空间中的甲烷分子增多，同时其扩散能力受到限制(压力越大，扩散系数越小)，且赋存的甲烷分子数量随着孔径的增大而增多；孔隙壁面附近存在明显密度峰值，峰值随着压力的升高而增大；采用外力法模拟施加压降，外力的增大会提升甲烷在孔隙中的流速，且赋存空间(孔径)越大，甲烷在孔隙中的流速越大。在开采过程中制造较大的压差，更有利于甲烷从页岩储层中脱离出来。

Abstract: Adsorption of gas is a primary method of shale gas storage in shale reservoirs. Furthermore, adsorbed gas constitutes the predominant source of shale gas production in the late stage. The accurate assessment of shale gas reserves is contingent upon a comprehensive understanding of the characteristics of methane storage in shale reservoirs under conditions of high temperature and high pressure. A molecular model of the organic matter-quartz composite pore system composed of quartz, casein and composite pores was constructed on the basis of the inorganic geochemical characteristics of shale in the Wufeng Formation (Longyi 1 subsection of the Sichuan Basin). The micro-existence characteristics of methane molecules in the composite pore system, and their flow-diffusion characteristics, were investigated by means of a range of simulation methods, including giant canonical Monte Carlo (GCMC), molecular dynamics (MD) and nonequilibrium molecular dynamics (EFMD). The results demonstrate that the number of methane molecules in the pore space increases with an increase in pressure, while their diffusion ability is limited (the higher the pressure, the smaller the diffusion coefficient), and the number of methane molecules increases with an increase in pore diameter. It is evident that there is a conspicuous density peak in proximity to the pore wall, which exhibits an increase in magnitude with an escalation in pressure. The pressure drop is simulated by the method of applying an external force, and an augmentation in the external force will result in an enhancement of the flow rate of methane in the pore space. Concurrently, the flow rate of the stored space will be increased by the external force. The pressure drop is simulated using the external force method, and an increase in external force will enhance the flow rate of methane in the pore space. Furthermore, it has been demonstrated that a larger endowment space (pore diameter) will result in a larger flow rate of methane in the pore space. It has been demonstrated that creating a larger pressure drop during the extraction process is more favourable for methane to detach from the shale reservoir.

1. 引言

“十三五”以来，中国石化在页岩气领域取得了一系列重大进展，其中四川盆地涪陵页岩气田的成功发现和探索性实践，成为国内页岩油气高效勘探的典范[1]。随着时间推移，页岩气勘探开发的重心逐步由浅部转向深部，并在实践中获得了重要发现与技术突破。据研究，四川盆地及周缘页岩气地质资源量21.9 × 10¹² m³，其中埋深3500~4500 m的深层页岩气占51%，为11.3 × 10¹² m³ [2]。2024年，涪陵页岩气田阳春沟区块新增探明储量1213.56 × 10⁸ m³，涪陵气田累计探明储量超万亿立方米，持续引领我国页岩气大规模增储和高质量勘探[3]。

页岩气是指赋存在富含有机质的成熟暗色泥页岩或高碳泥页岩中的天然气，其来源可为生物成因、热解成因或二者混合成因。该类天然气因有机质吸附作用及岩石裂缝与基质孔隙的存在而得以聚集和保存，赋存形态主要包括游离态、吸附态和溶解态(溶解态主要存在于有机质与地层水中)。与中浅层(<3500 m)页岩气相比，深层(>3500 m)页岩气埋深更大，页岩储集性能、赋存状态和含气性与中浅层相比有较大的差别[4]。目前，等温吸附实验仍是评价页岩储层气体含量的重要室内手段之一。然而，其测试条件通常限制在压力不超过40 MPa、温度低于90℃，与深部储层的高温高压地质环境存在明显差异，计算结果务必会与实际储量存在较大误差。分子模拟方法在近几年被引入页岩气研究，并取得了良好应用效果，分子模拟方法可以从分子角度认识和解释气体的动用机理，表征不同赋存状态页岩气的动态变化过程，常用模拟方法有巨正则蒙特卡洛(GCMC)和分子动力学(MD)方法。卢双舫等[5]基于GCMC方法模拟了甲烷和二氧化碳在伊利石狭缝内的吸附，认为范德华力普遍存在气固分子之间，当含有极性气体分子时还存在库仑力。ONAWOLE等[6]利用密度泛函理论(DFT)和MD方法研究了甲烷在二氧化硅–高岭石表面上的吸附行为，揭示了甲烷吸附为物理吸附，且甲烷与高岭石间的吸附强于与二氧化硅间的吸附。YU等[7]通过GCMC方法和MD方法研究了干酪根上的吸附和由吸附引起的收缩机理，指出了体积应力与吸附量之间呈线性相关。马欣健等[8]通过GCMC方法和MD方法揭示了页岩储层中注烟道气置换CH4的微观机制。

当前针对页岩气赋存的研究大多采用单一矿物构建分子模型，如干酪根分子模型、黏土矿物分子模型和无机矿物模型，但从常规扫描电镜实验可知，页岩储层中的孔隙较为复杂(粘土矿物粒间孔、脆性矿物粒间孔以及有机质孔等)，且这些孔隙是由有机质与黏土矿物和非黏土矿物之间相互粘结填充形成。针对这一问题，本文以四川盆地五峰组——龙一₁亚段页岩为研究对象，利用X-射线衍射实验与分子模拟相结合的手段，建立了有机质——石英复合孔隙模型，利用GCMC、MD和EFMD方法研究甲烷分子在该复合孔隙体系中的微观赋存特征及其扩散流动特征。望该项研究帮助深入理解页岩气在四川盆地页岩储层中的赋存状态与扩散流动特征。

2. 模型构建及计算方法

2.1. 模型构建

在构建分子模拟模型之前，对所取四川盆地五峰组——龙一₁亚段页岩样品(平均深度4174 m)进行全岩X-射线衍射分析，实验测试结果如图1所示。根据分析测试结果可知，矿物含量中脆性矿物含量高，其中石英的占比最高，含量范围在59.5~75.5%，平均值为68.7%；黏土矿物成分占比位居第二，含量范围在3.2~12.1%，平均值为9.4%。由此可知，四川盆地五峰组——龙一₁亚段页岩脆性较高。

Figure 1. Whole-rock quantitative analysis of shale samples from the Sichuan Basin

图1. 四川盆地页岩样品全岩定量分析

本文采用MS和Lammps软件协同建模与计算。页岩有机质的主体骨架是干酪根，本文基于Ungerer等人[9]提出的II型干酪根(C₂₅₂H₂₉₄O₂₄N₆S₃，图2(a))开展模拟研究，干酪根模型构建具体过程如下：(1) 对二维结构进行几何优化。选择COMPASS力场，使用Smart算法以Fine的计算精度运行几何优化(优化过程均选择Fine精度)；(2) 对几何优化后的二维结构进行退火模拟，控温函数选择Nose算法，总步数100,000步；(3) 进行分子动力学(MD)计算，选择NVT系综计算，步长1 fs，模拟时间设置为50 ps，得到能力最低的二维结构；(4) 将10个二维单元结构通过AC模块创建周期性结构，为了避免在单元结构的重叠，周期性结构的初始密度设置为0.4 g/cm³；(5) 对得到的初始结构进行NVT模拟，控温函数选择Anderen算法，随后逐级降温进行系综NPT模拟，控压函数选择Berendsen，步长1fs，模拟时间设置为100 ps。经过上述优化与驰豫后，得到密度为1.119 g/cm³的干酪根模型(图2(b))。由于干酪根的类型与成熟度的不同，真实的页岩中基质密度在1.1~1.40 g/cm³ [10]，由此可得构建的干酪根模型是合理的，可用于后续复合孔隙模型的建立。

(a) (b)

(c)

Figure 2. Kerogen unit and kerogen-group models. (a) Kerogen unit model; (b) Kerogen-group model; (c) Kerogen density variation curve

图2. 干酪根单元和干酪根基团模型。(a) 干酪根单元模型；(b) 干酪根基团模型；(c) 干酪根密度变化曲线

由X-射线衍射实验可知研究区页岩石英含量较高。对此参考美国矿物学家晶体结构数据库[11]导入石英晶体模型。为了尽可能与实际情况相符合，石英片层由模拟单元8 × 6 × 1扩胞形成，将石英片层模型和干酪根基质块体模型沿着Z方向组合，最终得到由干酪根和石英组成的有机质——石英复合孔隙模型，如图3所示。

Figure 3. Pore model of organic matter—non-clay mineral composite

图3. 有机质——非黏土矿物复合孔隙模型

2.2. 合理性验证

分子模拟结果依赖于所采用的分子结构与力场参数，因此其可靠性需通过实验加以验证，图4为实验结果与模拟结果对比图，分子模拟与实验结果尽管有一定差别，但处于相同数量级，且误差控制在10%以内，可以说吻合较好。造成误差因素可能是由于：在体积法等温吸附实验中，自由体积由氦气标定，然而He分子的动力学直径小于CH₄分子；而实验中自然矿物样品中可能存在杂质，其结构也并非理想结构，在类质同象替换、离子充填和表面晶格缺陷方面均可能与模拟单元之间存在差异，进而造成二者吸附能力的差异[4]。

Figure 4. Comparison of experimental and simulated excess adsorption capacities

图4. 实验与模拟超额吸附量对比图

2.3. 模型方法

2.3.1. 赋存模拟

等温吸附实验测试是测定每一个平衡压力点下的吸附量，Sorption模块对有机质–石英复合孔隙体系开展甲烷等温吸附模拟。模拟计算基于COMPASS力场，静电作用和范德瓦尔斯作用力分别使用Ewald求和方法与原子式计算方法。针对每一测试压力点，先进行1 × 10⁶步吸附平衡运算，再通过1 × 10⁷步采集平衡阶段的吸附量数据。复合孔隙吸附稳定构型可以获得甲烷分子在狭缝孔中瞬间构型图；通过统计计算求平均得到甲烷分子在孔隙体系中的局部密度分布。

2.3.2. 流动模拟

体系在x、y、z三个方向均设置周期性边界条件，先通过平衡分子动力学实现弛豫。当体系稳定后，依据爱因斯坦扩散理论，通过均方位移(MSD)计算扩散系数D，用以表征流体分子在纳米孔隙中的扩散性能。其计算公式为

$D=\frac{1}{6N_{\alpha }}\underset{t\to \infty }{\lim }\frac{\text{d}}{\text{d}t}{\displaystyle \sum _{i=1}^{N_{\alpha }}\langle {\left[r_i\left(t\right)-r_i\left(0\right)\right]}^2\rangle }$ (1)

系统对外部扰动的响应通常是非平衡热力学过程。因此，伊利石纳米狭缝中压力驱动的甲烷解吸采用外力方法处理，这是分子动力学模拟中常用的方法之一[12]。在该方法中，沿着流动方向(x)对每个甲烷分子施加一个额外的恒定力(F)，以模拟压力驱动条件。对于给定的压力差ΔP，相应的F通过以下公式计算

$F=\frac{1}{n}\frac{\partial P}{\partial L}=\frac{\Delta P}{n{L}_x}=\frac{\Delta P\cdot L_p\cdot H}{N}$ (2)

在压降的驱动下，分子将通过狭缝移动，其速度不仅取决于压力降，还取决于它们与伊利石表面的距离。在哈根–泊肃叶方程中[13]，不可压缩粘性流体的速度由下式给出

$\nu \left(z\right)=\frac{\partial P}{\partial L_x}\cdot \frac{1}{2\mu }\cdot z\left(H-z\right)$ (3)

其中$\frac{\partial P}{\partial L_x}$ 为压力梯度，$\mu$ 、H、z分别为流体的黏度、孔隙直径以及流体与壁面之间的距离，在箱式方法中，沿z方向的狭缝空间被划分为厚度为0.05 nm的箱。箱内甲烷分子的平均速度[14] [15]估算得出

$v_x=\frac{{\displaystyle \sum _{i\in \text{bin}}{m}_i}{v}_{i,x}}{{\displaystyle \sum _{i\in \text{bin}}{m}_i}}$ (4)

其中V_i,x是平板中甲烷分子的速度，m_i是其质量。

3. 结果和讨论

3.1. 赋存状态

为研究深层页岩储层高温高压条件下甲烷的赋存状态，本文参考王强等人川东南下志留统龙马溪组深层页岩气藏地层温度和压力参数[16]，设置模拟温度为110℃、压力为70 MPa，进行了甲烷在复合孔隙中吸附与扩散的模拟研究。并设置压力为5 MPa、30 MPa和70 MPa；孔径为3 nm、4 nm和5 nm的模拟情况，进行对比分析。

图5展示了不同孔径模型在不同压力条件下甲烷分子的赋存状态。从复合模型吸附构型图可知，随着压力的升高，吸附在孔隙壁面和游离在孔隙间的甲烷分子数量在显著增加。且孔径的增大会使得甲烷的储集空间增大，这也就使得在同一压力情况下，孔径大的孔隙模型中甲烷分子数量更多。

(a) H = 3 nm、P = 5 MPa (b) H = 3 nm、P = 30 MPa

(c) H = 3 nm、P = 70 MPa

(d) H = 4 nm、P = 5 MPa (e) H = 4 nm、P = 30 MPa

(f) H = 4 nm、P = 70 MPa

(g) H = 5 nm、P = 5 MPa (h) H = 5 nm、P = 30 MPa

(i) H = 5 nm、P = 70 MPa

Figure 5. Methane occurrence states in composite pores under different pore sizes and pressures

图5. 复合孔隙中不同孔径、压力条件下甲烷的赋存状态

孔隙内气体分子的分布特征是揭示页岩气吸附行为的重要指标，分子模拟可直观观察吸附平衡后甲烷在纳米孔隙中的微观分布。如图6所示，甲烷密度分布曲线存在两个峰值，该峰值出现在复合孔隙模型的壁面位置，甲烷分子受到孔隙壁的吸附作用，使其在紧邻孔隙壁的位置形成强密度带，及强吸附层(吸附相)；孔壁外部两侧的微弱峰值，是甲烷分子填充于干酪根内部的微孔中，将其称之为溶解相甲烷，如图7所示；吸附层之外，形成密度值小于吸附相的弱吸附层，其原因在于甲烷分子与孔隙表面的距离增加时，流–固作用力显著减弱，导致通道中心区域分子受壁面效应影响较小，从而表现出稳定的密度分布。且随着孔径的增大，密度曲线峰值出现一定数量的增大，但孔隙中间的甲烷游离相密度基本不变，孔径的增大主要增大了孔隙内游离相的体积[17]。

(a) 3 nm复合孔隙通道内甲烷密度分布 (b) 4 nm复合孔隙通道内甲烷密度分布

(c) 5 nm复合孔隙通道内甲烷密度分布

Figure 6. Methane density distribution curves in composite pores under varying pore sizes and pressures

图6. 复合孔隙中不同孔径、压力条件下甲烷密度分布曲线

Figure 7. Schematic of methane molecules stored within kerogen

图7. 赋存于干酪根内的甲烷分子示意图

3.2. 扩散特性

孔隙中的游离气与吸附在孔隙表面的吸附气共同参与吸附渗流，期间在某一温压下存在着吸附气体和游离气体的相互转化，并非吸附态气体分子一直为吸附态、游离态分子一直为游离态[18]。孔壁对于气体的吸附能力大小表现为一种统计趋势。气体的传输机理可以分为自由气相滑脱流动、自由气相努森扩散、吸附气相表面扩散、吸附气相对自由气相的诱导传输[19]，甲烷分子在孔隙内的扩散能力可以通过扩散系数(均方位移曲线斜率)来表征，扩散系数越大说明甲烷分子的自由程度越高。

如图8，随着压力的升高甲烷分子在有机质——非黏土矿物复合纳米狭缝中的流动特性受到限制，致使其扩散系数在高压力条件下降低。在低压条件下(5 MPa)，4 nm孔径更接近努森扩散的理想条件，气体分子与孔壁的碰撞频率较高，且滑脱流动效应显著，分子与孔壁的相互作用更有利于扩散，导致4 nm孔径的扩散系数高于5 nm孔径[20] [21]。在30 MPa的高压条件下，气体的扩散机制从努森扩散逐渐向过渡流和普通扩散转变，在这种情况下，扩散系数不再简单地与孔径成正比[22] [23]。而在70 MPa高压下，

(a)

(b)

Figure 8. Methane diffusion characteristics. (a) Mean-square displacement (MSD) of methane; (b) Diffusion coefficient of methane molecules

图8. 甲烷分子扩散特征曲线/图。(a) 甲烷的均方位移；(b) 甲烷分子扩散系数

扩散机制逐渐从努森扩散向普通扩散过渡，扩散系数与孔径的关系变得复杂。在70 MPa的高压下，随着孔径的增大，扩散系数可能会增大，这主要是由于孔径增大提供了更多的自由体积，降低了扩散阻力，同时分子间的碰撞对扩散的抑制作用相对较小[24]。这种现象表明，气体在纳米孔隙中的扩散行为受到压力、孔径和扩散机制的共同影响，且在不同压力条件下表现出不同趋势。

3.3. 流动速度

页岩气藏吸附渗流过程中孔壁对于甲烷等气体分子具有较强相互作用从而产生边界滑移现象，如图9所示。一般来说，流速随ΔP的增加而增大。当ΔP = 0时，分子将保持随机热运动，在孔隙方向上的统计速度为零。当施加定向压降时，分子受到压力驱动，其速度沿该方向增加。图10展示了甲烷分子的速度随压降ΔP的变化情况。

Figure 9. Schematic of fluid flow in a nano-slit, velocity profile and the corresponding flow regions with slip length [25]

图9. 纳米狭缝中的流体流动以及速度分布及其滑移长度的相应流动区域示意图[25]

Figure 10. Variation of methane molecular flow velocity with pressure drop in a 3 nm pore

图10. 3 nm孔径中甲烷分子的流速随压降的变化关系图

页岩气藏吸附渗流过程中孔壁对于甲烷等气体分子具有较强相互作用从而产生边界滑移现象，如图9所示。一般来说，流速随ΔP的增加而增大。当ΔP = 0时，分子将保持随机热运动，在孔隙方向上的统计速度为零。当施加定向压降时，分子受到压力驱动，其速度沿该方向增加。图10展示了甲烷分子的速度随压降ΔP的变化情况。

在低压降(ΔP)时，分子受热运动主导，呈现出微弱的定向流动。当压降ΔP足够大以克服热运动时，定向流动占主导，其速度随压降ΔP增加。当施加的外力从1 MPa增加至40 MPa时，孔隙中间甲烷分子流速从2.39 m/s增至120.34 m/s，文献[26]-[28]也报道了因压力梯度增加而导致的类似速度增量。

图11展示了在固定压力(40 MPa)下，不同孔径狭缝中甲烷的速度分布。较宽狭缝中的分子往往具有更快的速度。对于H = 3 nm，速度较低且变化幅度较小，这是因为大部分甲烷分子吸附在孔隙表面，孔径越小受到表面的限制越明显。当H增加时，气相中的分子数量增加。在距表面较远的距离处，这些分子受到的表面吸引力较弱，在压力驱动下具有更快的定向运动，这表明孔径对窄纳米狭缝中甲烷流速有很大影响。其他作者也探讨了气体流动中的孔径效应，Yu等人[29]也发现相关的研究规律。

Figure 11. Methane molecular flow velocity curves in different pore sizes at a pressure drop of 40 MPa

图11. 压降40 MPa时不同孔径中甲烷分子的流速曲线图

4. 结论

首先本文采用分子动力学模拟技术构建了有机质——石英复合孔隙模型，通过等温吸附实验和模拟结果相比对保证所构建模型的可靠性。其次，探究了甲烷分子在不同深层温压条件下的赋存形态。最后，通过施加压降的方法使得赋存甲烷分子扩散。得到以下几点认识：

(1) 在相同孔径条件下，深层页岩因处于高温高压环境，其孔隙内甲烷分子含量显著高于中浅层页岩。甲烷密度剖面呈近似轴对称分布，孔壁附近存在明显峰值，且该峰值随压力升高而增强；而游离相密度在不同压力下变化不明显。

(2) 甲烷在纳米孔隙中的扩散系数随压力增加逐渐减小，高压情况下的扩散系数比低压情况下小的多；在低压条件下(5 MPa)，4 nm孔径滑脱流动效应显著，致使4 nm孔径的扩散系数高于5 nm孔径；高压条件下，压力升高导致孔隙内分子数增多，分子间碰撞频率加剧，从而抑制了扩散。

(3) 采用外力法模拟施加压降。在给定外力条件下，随着孔径增大，甲烷分子的平均流速增大。对于给定的狭缝，甲烷分子的流速随外力的增加而增大。孔径越大施加的外力越大，更有利于甲烷分子的解吸。

基金项目

重庆科技大学研究生创新项目(YKJCX2420135)。