多孔介质内甲烷碳同位素迁移扩散的分子动力学模拟研究

刘 畅，张慧，赵耀中，陈威，2

(1.苏州大学能源学院，江苏 苏州 215006；2.中国航天科技集团有限公司氢能工程技术研发中心，北京 100074；3.北京航天试验技术研究所，北京 100074)

随着常规油气资源的日趋紧缺和出于能源清洁的考虑，页岩气这一清洁能源已成为多国学者研究的热点之一[1]。页岩气作为一种非常规能源，由于其具有储量丰富，燃烧时污染程度低的优良性能[2]，在能源领域发挥着越来越重要的作用。近年来，水平钻井和压裂技术的发展提高了页岩气的开采能力，促进了美国“页岩气革命”的成功。学者们对致密页岩内的气体流动展开了相关研究，但是目前对于页岩气在纳米级孔尤其是在含孔喉纳米孔中的吸附、扩散和同位素分馏特征仍不十分清晰[3]。

稳定碳同位素气体的组成、分馏可以为各成熟阶段的天然气提供有价值的地球化学信息，如气体的成因、示踪运移、聚集和散失等[4-7]，也被作为探索页岩内部气体分子微观尺度非线性输运过程的新方法。甲烷气体产出具有明显的多尺度流动的特征，包含了渗流、吸附/脱附和努森扩散等大量分子层次的非线性流动，其流体的渗透率和有效扩散系数等参数不再是常数[8]。在纳米级孔内不但存在烷烃气体分子的相对运动，而且存在甲烷同位素气体(12CH4和13CH4)的分馏(即在流动、吸附和热解等物理或者化学过程中，轻、重同位素分子的丰度在物相之间重新分配，12CH4和13CH4相对比例发生浮动)。由于实验技术和测量手段的限制，目前还无法直接观测热采工况下致密页岩内部甲烷的赋存状态和微观流动现象。因此，模拟成为研究同位素微观状态特征和扩散形式的有效手段。ZHANG等人使用分子动力学模拟的方法研究了页岩气在孔喉处的吸附特性，并且指出了页岩气藏中孔喉的下限[9]；NADER等人利用晶格玻尔兹曼模型研究了喉-孔比、Kn数对页岩气流动特性和滑移速率的影响[10]；一般认为甲烷碳同位素在受限空间中的分子扩散效应和物理吸附/脱附是造成其分馏的两个主要原因。XIA等人通过建立多孔流动模型研究吸附/脱附过程对同位素分馏造成的影响，LIU等人通过建立数学模型比较了烷烃碳同位素的分馏程度和吸附气比例的关系[11]。研究页岩气甲烷碳同位素的吸附和扩散过程能够强化对其流动特性的理解，有助获得更多页岩气的流动和储存信息；探究页岩气甲烷碳同位素在页岩气藏中的分馏程度对追踪页岩气的生成和演化、评价页岩气的相对吸附能力具有重要的意义。

因此，本文借助分子动力学的方法研究了甲烷同位素气体在页岩纳米级孔中(包含孔喉)的吸/脱附和扩散现象，获得不同工况下12CH4和13CH4扩散能力和分馏特征，旨在分析不同孔结构(包括孔喉大小)对甲烷同位素分馏的影响规律。

1 分子动力学理论基础与模型

1.1 分子动力学理论

本文应用分子动力学模拟软件(LAMMPS)研究甲烷碳同位素在有无孔喉存在的纳米级有机孔中的扩散速率和分馏特征。模拟甲烷碳同位素在压差的驱动下沿着碳纳米管流动。粒子间相互作用采用截断及移位的Lennard-Jones(LJ)势能函数[12]，即：

(1)

公式(1)中，E为分子势能，r为分子间距，rc为截断半径，本文中的rc=10 Å，σ是相互作用的势能正好为零的两体距离，ε是势能阱的深度。相同粒子间的相互作用势参数σ和ε[12]见表1，不同粒子间的σ和ε参数由Lorentz-Berthelot混合规则确定[13]：

表1 相同粒子间的相互作用势参数

(2)

1.2 同位素分馏模型

分馏模型采用两种不同结构的碳纳米管，用来模拟甲烷碳同位素的分馏。如图1(a)是含有孔喉存在的碳纳米管，其长度为20.4 nm，直径为4 nm，孔喉的长度为0.4 nm，直径为1 nm，孔喉左右两侧碳纳米管长度相同，均为10 nm。为了对比孔喉的影响，构建相同尺寸的不含有孔喉的碳纳米管用来研究甲烷同位素的分馏，即：碳纳米管的长度为20.4 nm，直径为4 nm(图1(b))。

图1 有/无孔喉存在的碳纳米管模型

为了比较甲烷同位素在不同孔内的分馏特征，12CH4和13CH4同位素在两种管内的扩散过程被分别模拟。首先使用巨正则蒙特卡罗方法，模拟了页岩在地下贮存353 K左右的温度和20 MPa的压力条件，使得模拟结果更加接近实际情况。在温度为353 K、20 MPa的压强下在含有孔喉的碳纳米管中生成一定数目的12CH4分子。模拟初始阶段，先通过Nose-Hoover控温器使管内气体分子在NVT系综下维持温度为353 K，经过时长为2 ns的弛豫时间后粒子温度、压力等热力学量到达稳定，且在一定时间内不发生波动，则系统达到热平衡，然后给每个粒子在流动方向施加外力，达到模拟孔隙内甲烷气体压差作用下流动的目的；13CH4粒子在同样的条件下生成并在相同压差条件下流动，分别统计并计算13CH4和12CH4的迁移扩散系数进行分析。本文考虑了非理想气体效应，利用Peng-Robinson状态方程求出逸度系数来表征实际气体偏离理想气体的程度[14]。在不含孔喉的情况下，以同样的方法生成粒子并附加相同的条件，分别统计并计算13CH4和12CH4的迁移扩散系数，与含有孔喉的情况进行对比分析，旨在探究孔喉的存在对甲烷碳同位素分子分馏特性的影响。

2 气体流动特性和扩散系数的计算

2.1 气体流动特性

与传统宏观的连续流体流动状态不同，由于高努森数效应的存在，达西方程已经不适用于描述甲烷同位素气体在纳米孔中的流动。页岩气的流动状态可由努森数来判定。努森数(Kn)的定义为[15]：

(3)

(4)

(5)

公式(3)中，λ为平均分子自由程，H为流动通道的特征参数。公式(4)中，R为气体常数，T为绝对温度，Mw为气体分子的摩尔质量，μ为粘度系数，P为气体压强。公式(5)中，KB为玻尔兹曼常数，dm为分子直径。当Kn<10-3时，流动属于连续流；当10-310时，流动形式为分子自由扩散[16]。本文模拟的Kn数的范围为0.06～0.18，属于滑脱流和过渡流的范围，与传统无滑移的连续流不同，需要考虑壁面的滑移效应。

2.2 扩散系数

按流动方向(z轴)统计区域均分成长度为4 Å的小区域，每个小区域通过的摩尔流率计算公式如下：

(6)

公式(6)中，A为碳纳米管的截面积，L为小区域的长度(此处为4 Å)，Ni为小区域内甲烷分子的个数。当系统达到稳定之后，由于边缘效应的存在，进口区域和出口区域的摩尔流率波动很大，但中间区域摩尔流率较稳定，因此选取稳定区域的摩尔流率平均值作为页岩有机纳米孔的平均摩尔流率，从而可求扩散系数[17]：

(7)

公式(7)中，D为12CH4的扩散系数，D′为13CH4的扩散系数。其中R为气体常数，T为温度，Lpore为碳纳米管的长度，F和F′分别为12CH4和13CH4分子的摩尔流率，ε为孔隙度，A为碳纳米管的截面积，ΔP是通过碳纳米管的压差，压差与外力的关系由下式给出：

(8)

公式(8)中，n为外力驱动的粒子数目，Fex为外加驱动力。由于本文中的驱动力以压差的形式表示，孔喉两端施加的10 MPa压力差为粒子扩散的驱动力。因此，12CH4和13CH4的迁移扩散系数的公式可以表示如下：

(9)

3 结果和讨论

3.1 甲烷同位素分子在孔内的径向密度分布

为了直观地反映纳米级孔内甲烷碳同位素分子的吸附特性和密度分布，图2给出了甲烷同位素分子在无孔喉存在的纳米通道的迁移扩散现象，图3给出了甲烷同位素分子在无孔喉模型中径向分布函数收敛后的相关的数据(压差为10 MPa)。甲烷同位素分子在靠近壁面处形成第一和第二吸附层，局部密度远高于孔中间分子游离区的密度。由于壁面与甲烷碳同位素分子之间存在较强范德华力的作用，粒子易在吸附位点取得平衡并形成吸附层。并且由于甲烷碳同位素分子数量多，在壁面处形成竞争吸附，由此形成第一吸附层和第二吸附层，靠近壁面的第一吸附层甲烷碳同位素分子的密度远大于第二吸附层甲烷碳同位素分子的密度。而当甲烷碳同位素远离壁面时，壁面对分子的吸附作用减弱，而分子间的作用力占主导。

图2 无孔喉模型中CH4在10 MPa压差下沿着Z轴方向的流动

图3 无孔喉模型中12CH4和13CH4在10 MPa压差下的径向密度分布

图4为甲烷同位素在含有孔喉的纳米孔中的流动模型。由图4可知，壁面处13CH4分子的密度更高，这是由于12CH4和13CH4在相同的温度条件下，质量大的13CH4粒子处于更深的势阱，更容易被壁面所束缚。

图4 甲烷同位素在含有孔喉的纳米孔中的流动模型

孔喉的存在通常会引起甲烷同位素密度分布在其前后产生差异。甲烷同位素在10 MPa的压差下沿着Z轴正方向迁移扩散。图5为甲烷通过孔喉前后的径向密度分布。从图5中可以看出，流过孔喉之前甲烷同位素的密度分布明显高于流过孔喉之后的密度分布，这是由于孔喉对甲烷碳同位素的通过具有一定的选择性，只有带有足够能量的分子才能越过孔喉的能量壁垒，进入到孔喉的另一侧[18]。而且由于孔喉处直径较小，一部分粒子会与孔喉附近壁面相撞而改变运动路径，从而降低了自身通过孔喉的概率。

图5 甲烷通过孔喉前后的径向密度分布

3.2 甲烷同位素分子在孔内的速度分布

在本文模拟中，甲烷碳同位素分子(12CH4和13CH4)在纳米通道中的流动形式处于滑脱流和过渡流两种流态(0.06

此外，对比图6(a)、图6(b)可知，甲烷碳同位素分子在含孔喉有机孔中流动的径向分布速度小于无孔喉有机孔中的速度。这是由于气体分子会与流动方向上的壁面发生碰撞，粒子存在一定的概率不能通过孔喉。而且由于孔喉部分直径小，粒子可以自由移动的区域减小，由此阻滞了甲烷碳同位素分子的迁移扩散过程，从而降低了径向分布速度。

图6 12CH4和13CH4分子在有无孔喉存在的碳纳米管内流动的径向速度分布

3.3 同位素迁移扩散系数和分馏程度

页岩气随着生产过程的进行，孔内压力逐渐降低，由于甲烷同位素分子脱附能力的差异导致其分馏程度发生变化[19]。通过改变纳米孔中对甲烷的加力大小，可以模拟甲烷同位素在不同压差下的迁移扩散。图7(a)、7(b)分别表示在有孔喉和无孔喉页岩有机孔中的甲烷碳同位素分子迁移扩散系数随压差的变化。从图中可以看出，随着压差降低，12CH4和13CH4的迁移扩散系数均减小，并且有孔喉存在时，由于孔喉的阻碍和内部的吸附作用，甲烷碳同位素分子的迁移扩散系数更小，即迁移扩散速度减慢。此外，在任一压差下12CH4的迁移扩散系数均大于13CH4，这是由于13CH4分子的质量大，碳纳米管对其吸附能力更强；此外，受驱动力影响，质量大的13CH4的扩散速度更慢，由此造成同位素分馏。

图7 甲烷碳同位素分子迁移扩散系数随压差的变化

吸附和扩散是导致同位素分馏两个主要原因，本文采用13CH4和12CH4扩散系数的比值δ表征同位素分馏程度的大小。图8为甲烷碳同位素分子的分馏程度随压差的变化。由图8可知，在无孔喉有机纳米孔中，12CH4和13CH4存在分馏现象，并且分馏程度小于同位素在含孔喉有机孔中的分馏程度。本文中孔喉直径(1 nm)与甲烷分子直径(0.38 nm)比较接近，孔喉对甲烷分子的吸附作用明显较大，两种同位素分子可以自由移动的区域较小。因此当两种粒子进入孔喉时，受范德华力的影响较大，并且由于13CH4的质量大于12CH4，所处的势阱位置更深，更易吸附在壁面，所以与无孔喉有机孔相比，在孔喉纳米孔中甲烷碳同位素分子的分馏程度更大。

图8 甲烷碳同位素分子的分馏程度随压差的变化

甲烷碳同位素气体在有无孔喉的纳米孔中的分馏程度均随压差的降低先增大后减小，即甲烷碳同位素分子先变 “轻” 后变 “重”，模拟结果与LI等人获得的实验结果变化趋势一致[20]。这是因为随着压差降低，两者中吸附能力较弱的12CH4首先从壁面大量脱附，而13CH4的脱附程度较弱，由此造成两种同位素分子分馏程度增大。继续降低压差，吸附能力更强的13CH4在压差的作用下开始从壁面大量脱附，因此导致两种分子的分馏程度逐渐减弱[21]。

4 结语

本文应用分子动力学模拟方法研究了353 K和不同压差(10～30 MPa)下12CH4和13CH4分别在包含和不包含孔喉的页岩有机纳米孔中的径向密度分布以及摩尔流率，分析了孔喉结构对甲烷碳同位素(12CH4和13CH4)分馏程度的影响，得出以下结论：

(1)当有机孔不含孔喉时，甲烷在孔壁面的吸附和扩散作用会引起甲烷碳同位素分馏。在相同压差下，12CH4的质量更小，因此12CH4有更高的扩散速度，同时壁面对13CH4更强的吸附能力进一步减弱13CH4的扩散速度，导致了碳同位素气体的分馏。

(2)当有机孔含有孔喉时，由于孔喉的阻碍作用，气体分子会与流动方向上的壁面发生碰撞，降低通过孔喉的概率，从而导致有孔喉存在时同位素的径向分布速度变小。

(3)由于孔喉直径小，使得含孔喉的纳米孔中绝大部分区域属于吸附区，因此进入孔喉的粒子受壁面吸附影响，而13CH4的质量大于12CH4的质量，所处的势阱更深，导致孔喉对13CH4的截留作用更强，进而增加了同位素分馏程度。

(4)由于同位素吸附/脱附和扩散能力的差异，随着压差的降低，页岩气甲烷碳同位素分馏程度在不同结构的纳米孔内均为先增大后减小。