混合无机纳米孔内流体运移规律研究*

2023-12-22高明玉刘丹董志立梁爽王京坤罗鑫

能源化工 2023年5期

高明玉，刘丹，董志立，梁爽**，王京坤，罗鑫

(1. 东北石油大学石油工程学院，黑龙江大庆 163318；2. 东北石油大学提高油气采收率教育部重点实验室，黑龙江大庆 163318；3. 国家管网集团油气调控中心，北京 100013；4. 大庆油田第八采油厂，黑龙江大庆 163514；5. 大庆油田第四采油厂，黑龙江大庆 163511)

随着社会经济发展，对油气资源的需求日益增大，非常规油气资源越来越受到重视，页岩油也随之成为世界各国油气勘探的焦点。早在2012 年，美国便依靠大力开发页岩油实现了原油产量剧增，2019 年美国超过俄罗斯和沙特阿拉伯成为全球最大的原油生产国。由此可见，页岩油的勘探开发对世界能源格局有着不可忽视的影响[1]。

页岩油是以页岩为主的页岩层系中所含的原地滞留油气资源，圈闭界限不明显，无法形成自然工业产能[2]。与常规储层不同，页岩储层具有极低的渗透率和丰富的纳米孔隙，对于几纳米到几十纳米的微小孔隙，传统流体力学理论和试验方法难以描述流体的吸附、运移等细节，基于牛顿力学的分子动力学模拟则能很好地解决这些问题[3]。分子动力学模拟以分子为研究对象，将整个系统作为有一定特征的分子集合，利用牛顿运动定律对系统内分子的运动轨迹进行模拟，并通过系综计算体系的结构和性质。尽管页岩油气赋存及运移规律的研究已经取得了诸多成果，但页岩储层矿物成分复杂，存在大量多种成分混合的纳米孔，对于混合孔隙内流体运移规律的相关研究较少，仍需进一步分析。

笔者以辛烷和水为研究对象，采用分子动力学模拟方法开展了石英和伊利石混合的无机质孔隙内流体赋存状态和运移规律的研究，通过分析孔隙内流体的密度分布、速度分布、滑移长度、黏度等参数，分别讨论了孔径大小、温度、压力梯度及含水率对纳米孔隙内流体赋存和运移的影响。

1 模型构建及模拟方案

1.1 模型构建

页岩基质骨架由有机质和无机质组成，有机质主要成分为干酪根，无机质主要成分为石英和黏土矿物，黏土矿物中以伊利石、绿泥石及伊蒙混层为主要成分，部分含有极少量绿蒙混层，其中伊利石占比最大，平均含量达到74.1%[4]，因此，以石英与伊利石组合模型构建页岩孔隙。伊利石是以Al原子为中心的八面体层，夹在2 个以Si 原子为中心的四面体之间作为2∶1 型黏土矿物。在每个复制单元晶胞中，一个Si4+被一个Al3+取代，形成同构取代。层间反阳离子(钾离子，K+)在伊利石层间空间中随机分布，以抵消伊利石层内同构取代引起的静电电荷，伊利石结构原子坐标及晶胞参数分别见表1、表2。石英是一种由SiO2四面体通过共同顶点连接成的三维骨架结构。在石英晶体中，硅原子位于四面体中心，1 个硅原子与4 个氧原子成键，同时1 个氧原子与2 个硅原子成键，氧原子处于硅氧四面体的公用顶点上。

表1 伊利石模型原子坐标

表2 伊利石模型晶胞参数

晶体的形貌是晶体晶面生长发育的结果，在晶体生长发育的过程中，生长速度快的晶面往往消失，生长速度慢的晶面容易保留下来，成为晶体形貌最主要的晶面。在不破坏Si—O 键的情况下，伊利石的(001)面是伊利石形貌最重要的晶面，该晶面用于考察外界与晶面的相互作用，有足够的代表性。Chilukoti 等[5]通过试验发现液态烷烃更倾向于吸附在α-SiO2的(100)晶面。因此，模型构建将α-SiO2的(100)晶面作为SiO2的代表面，与伊利石的(001)面组合成伊利石-二氧化硅模型，模型结构见图1；狭缝孔隙由两个平行的伊利石-二氧化硅组合壁面构成，向两个壁面中添加一定数量的n-C8H18[6-7]，构成孔径分别为2，5，8 nm 宽的孔隙。模型在3个方向均采用周期性边界条件，为了消除边界效应，在3 个方向上分别添加了0.1 nm 的真空区域，伊利石-二氧化硅混合页岩孔隙内流体分布初始结构见图2。

图1 伊利石-二氧化硅组合模型结构

图2 伊利石-二氧化硅混合页岩孔隙内流体分布初始结构

1.2 模拟方案

为了保证模拟结果的准确性，分子动力学模拟过程中选取的势能参数至关重要。采用SPC/E 模拟水分子，ClayFF 力场模拟石英和伊利石基质[7-8]，此力场广泛应用于黏土界面模拟。ClayFF力场仅考虑水分子、羟基、可溶性多原子分子和离子间的键伸缩和键角弯曲项，其他所有的相互作用力均采用非键相互作用势即Lennard-Jones势与库仑力之和进行模拟。孔隙内的n-C8H18分子采用OPLS-AA(Optimized Potentials for Liquid Simulation)全原子力场模拟[7]，其特点是对于不同原子间的势能参数通过试验结果拟合得到，此力场被广泛应用于烷烃、聚合物、生物大分子等，具有很高的可靠性。短程非键范德华相互作用使用Lennard-Jones(12-6)表示，其截断半径设置为1.2 nm，符合截断半径不能大于最小模拟盒子尺寸一半的规则。采用基于傅里叶的Ewald求和方法计算了长程静电力相互作用，粒子与粒子间的粒子网格采用PPPM方法，精度值设置为10-6。通过Parrinello-Rahman稳压器控制体系的压力，通过Nose Hoover恒温器控制温度，模拟初始在等压等温NPT系综下进行，温度360 K，压力30 MPa，体积波动仅设置在z方向。设置模拟的时间步长为1 fs，在NPT系综下弛豫1000 ps，取最后500 ps数据用来分析。

2 结果分析

2.1 孔喉尺寸对流体运移影响

温度为360 K，孔径分别为2，5，8 nm 狭缝内流体的连续密度分布曲线及压力梯度为0.002 kcal/(mol·Å)(换算单位为MPa/nm)的速度分布曲线见图3。

图3 不同孔径的密度(速度)分布曲线

由图3 可见：不同孔径内流体均为对称分布，孔径2 nm 狭缝内两侧各形成2 个吸附层，孔径分别为5 nm 和8 nm 狭缝内两侧均形成四个吸附层，且各吸附层峰值密度逐渐减小；吸附相密度与狭缝孔径大小有关，2 nm 孔径中第一吸附层密度约为0.81 g/cm3，5 nm 孔径和8 nm 孔径中第一吸附层密度相近，约为0.92 g/cm3。由此可见，随着孔径宽度的增大，吸附层数量、第一吸附层密度均呈先增大后趋于稳定的趋势。

值得注意的是，孔道边界处流体受到壁面吸附作用以及z 方向的驱动力，而孔道中央处流体仅受到驱动力作用。孔隙内流体分布及受力情况见图4。

图4 孔隙内流体分布及受力情况

由图4 可见：孔径分别为5 nm 和8 nm 狭缝内存在吸附相流体和体相流体，而2 nm 狭缝内充满吸附相流体，不存在体相流体，这是由于在较小的孔隙中，无机质表面对烷烃分子的作用力强于烷烃分子内部的作用力，因此所有烷烃均吸附在无机质表面；当孔径增大时，无机质表面对烷烃的引力逐渐减弱，因此在距离壁面较远的孔道中央形成体相流体，呈游离态，密度为0.65 g/cm3，与NIST 提供的试验数值一致。

根据图3 中不同孔径狭缝内流体的速度分布，其中散点为模拟结果，实线为该散点的抛物线拟合结果，孔径分别为2，5，8 nm 狭缝内流体流速峰值分别为0.00029，0.00178，0.00599 Å/fs，结果表明在其他条件相同时，与小孔径相比，大孔径内流体流速远大于小孔径内的流体流速。狭缝内流体并非以同一速度流动，受到无机质壁面的影响，壁面附近流体流速低，孔道中央流体流速高，速度分布为抛物线形状，可由二次函数拟合：

式中：v为流体流速，z为z方向x轴上距中心点的距离，a和c为拟合参数。

两光滑平面内不可压缩层状流体的稳态速度剖面v(z)可通过Poiseuille 方程描述：

式中：n表示分子数，F表示驱动力，η表示流体黏度，W表示孔径宽度。

根据连续流体力学理论假设，流体在固体壁面处速度为0，但烷烃在狭缝内流动时，固液界面存在滑移现象[8]，在纳米级别体系中无法忽略同处在纳米级别的滑移长度，因此引入滑移长度Ls建立更为通用的边界条件，狭缝内具有滑移的Poiseuille 流动的速度剖面为：

烷烃黏度可表达为：

滑移长度定义为从无机质壁面到流体速度为0 处的外推距离[9]，在流体力学中，边界条件模型包括正滑移模型、无滑移模型以及负滑移模型[10]，滑移长度模型见图5。

图5 滑移长度模型

图5 中，vs为滑移速度，b 为滑移长度。b ＞0 时为正滑移模型，表示狭缝内所有流体均参与流动，且壁面处流体流速不为0；b=0 时为无滑移模型，表示壁面流体流速恰好为0；b＜0时为负滑移模型，表示壁面附近具有厚度为b 的黏滞层，该部分流体不参与流动。

不同孔径的流体滑移长度及黏度见图6。

图6 不同孔径的流体滑移长度及黏度

由图6 可见：该驱动力条件下滑移长度均为负值，随孔径增大，负滑移长度先是急剧减小，后趋于平缓，这是由于小孔径内两侧壁面距离很近，吸附在一侧的流体仍会受到另一侧壁面的吸引从而减小了单侧壁面对流体的相互作用，随着孔径增大，受另一侧壁面的影响减小，不参与流动的黏滞层厚度也因此增大。

针对流体黏度，分别对比狭缝内全部流体平均黏度与体相流体黏度。ηeff为平均黏度，ηcenter仅对中间体相流体拟合，为体相流体黏度。2 nm 孔径内流体均为吸附相，不存在体相流体，因此ηeff和ηcenter基本一致。孔径分别为5 nm 和8 nm 狭缝中存在体相流体，由于烷烃的剪切稀释作用导致ηeff小于ηcenter。孔径为2，5，8 nm 狭缝内流体的平均黏度分别为0.24，0.30，0.33 mPa·s，说明随着孔径的增大，流体的平均黏度有增大的趋势，逐渐接近体相流体的黏度。

2.2 压力梯度对流体运移影响

驱动力分别为0.002，0.0025，0.003 kcal/(mol·Å)，孔径为5 nm 狭缝内流体的速度分布曲线见图7，不同压力梯度下流体的滑移长度及黏度见图8。

图7 不同压力梯度下流体的速度分布曲线

图8 不同压力梯度下流体的滑移长度及黏度

由图7 和图8 可见：不同驱替压力梯度下孔道内流体的密度分布基本重合，因此受限空间内烷烃的赋存特征与驱动力无关。随压力梯度的增大，流体流动速度增大，滑移长度也增大，压力梯度从0.002 kcal/(mol·Å)增大到0.0025 kcal/(mol·Å)时，滑移长度由-0.35 nm 增大到-0.13 nm，这是由于不流动的黏滞层中远离壁面的部分逐渐被驱动，导致黏滞层厚度降低。压力梯度从0.0025 kcal/(mol·Å)增大到0.003 kcal/(mol·Å)时，滑移长度由-0.13 nm 增大至0.004 nm，此时狭缝中不存在黏滞层，包括靠近壁面处的流体在内，狭缝内全部流体参与流动。在三种压力梯度下，流体黏度几乎不变，均在0.30 mPa·s 左右，可见狭缝内流体的平均黏度几乎不受驱动力的影响；不同矿物孔隙中流体表观黏度也有差异。

Zhang 等[11]研究了不同压力梯度下方解石孔隙中辛烷黏度约为0.48 mPa·s，王森等[6]在石英孔隙中研究了不同压力梯度下辛烷黏度约为0.29 mPa·s，与笔者的黏度结果十分接近，原因是混合无机孔隙中石英对烷烃的相互作用起到了较大的影响。

2.3 温度对流体运移影响

不同温度下第一吸附层峰值密度及体相流体密度见图9。

图9 第一吸附层峰值密度及体相流体密度

由图9 可见：随着温度的升高，第一吸附层峰值密度由0.92 g/cm3变为0.81 g/cm3，有明显降低趋势，即升温更有利于烷烃分子从无机质表面逃离，壁面对烷烃的吸附能力受到抑制，因此吸附相密度减小。不同温度下流体的密度(速度)分布曲线见图10，不同温度下流体的滑移长度及黏度见图11。

由图10 和图11 可见：升温过程中，流体流速增大，平均黏度降低，流动有所增强；滑移长度由-0.35 nm 增大至0.8 nm，这意味着360 K 温度下，狭缝两侧存在厚度为0.35 nm 的黏滞层；随着温度升高，滑移长度逐渐由负变为正，黏滞层厚度在逐渐降低，边界处原本不流动的黏滞层开始慢慢被驱动，当滑移长度变为正值时，狭缝内全部流体均参与流动，因此热采可成为提高页岩油采收率的有效措施[12]。

2.4 含水率对流体运移影响

含水率分别为20%，50%，80%的油水平衡构型见图12。

由图12 可以直观地看出当孔隙中存在油水两相时，石英-伊利石混合壁面优先吸附水，含水较低时，在两侧壁面附近形成水膜，孔道中央为烷烃分子，形成水-油-水的夹层结构，随着含水率的增大，壁面附近的水膜厚度增大，当含水率达80%时，混合无机纳米孔内出现水桥，烷烃分子在孔道中央聚集成团簇状。

不同的含水率条件下流体的密度及速度分布见图13～15。

图13 20%含水率时流体密度及速度分布

图14 50%含水率时流体密度及速度分布

图15 80%含水率时流体密度及速度分布

由图13～15 可见：含水率为20%时，水在壁面处仅形成一个吸附层，吸附层峰值密度为1.101 g/cm3；当含水率增加到50%时，吸附层数量增加至4 个，第一吸附层峰值密度增大至1.26 g/cm3；含水率继续增大至80%，吸附层数量和吸附层峰值密度趋于稳定，几乎不再变化。这是由于含水率很低时，壁面吸附的水层并未达到饱和，随着含水率的增加，吸附水层的数量、厚度也会随之增加，当吸附达到饱和时，吸附水层的数量、厚度则趋于稳定。观察烷烃的密度分布可以发现，含水率20%时，烷烃也存在微小的吸附效果，这是由于吸附水层厚度很小，不足以屏蔽壁面和烷烃间的相互作用，这种相互作用随吸附水层的增大而减小，含水率50%时这种相互作用被水层完全屏蔽[13]。在含水率不断增加的过程中，体相油区越来越小，油水两相区域逐渐增大，当含水率达80%时，不再存在体相油区，此时边界处为吸附水相，狭缝中其余空间均为油水混相区域。此时油水速度分布均呈抛物线形，Poiseuille模型拟合结果表明，油水速度分布呈明显的不连续性，即烷烃和水之间存在液-液滑移现象。在第一吸附层附近水流速为0，这代表了混合无机纳米孔表面的负滑移现象，存在一定厚度不参与流动的黏滞水层，滑移长度及油相黏度。不同含水率条件下烷烃滑移长度及油相黏度见表3。