气驱油油气混相过程的界面传质特性及其分子机制

俞宏伟，李实，李金龙，朱韶华，孙成珍,＊

1中国石油勘探开发研究院，提高石油采收率国家重点实验室，北京 100083

2中国石油吉林油田分公司，吉林 松原 138099

3西安交通大学，动力工程多相流国家重点实验室，西安 710049

1 引言

提高原油采收率一直是石油工程领域的一个热点话题。一次、二次采油仅能采集油藏中30%-40%的油，还有大量的原油难以脱附下来，推动采油技术得到了不断发展1-5。近年来，注气驱油(包括CO2、N2以及烃气)6-8在二次采油以及三次采油过程中发挥了重要的作用；尤其是二氧化碳驱油技术作为将提高采收率与CO2封存同时实现的途径，受到了越来越多的关注9-14。气驱油技术提高原油采收率的主要机理是驱替相气体可以很好地溶解于原油中14,15，原油膨胀导致界面张力、原油粘度与密度都降低，促进原油在油藏中的流动，进而提高石油的采收率。

气驱油过程涉及到多种物理化学问题，如油气两相物理性质的变化、界面特性的演变、孔隙内的油气多相流动等16-21。近年来针对于气体在油相中的溶解行为，Zhang等人22和Yang等人23对CO2在油中的分散特性进行了研究，前者从模拟角度出发对CO2溶解度与辛烷密度降低/溶胀之间的关系提供了定量的认识，后者利用实验手段对超临界条件下CO2-烷烃的体积膨胀规律进行了分析。也有学者通过分子动力学(MD)模拟研究超临界状态下的气体对水-油体系界面和输运特性的影响24。同时，部分学者发现外部环境的改变会导致油气两相状态发生变化，如Kiran等人25发现温度对原油的溶解度和粘度有重要影响，Cao等人16和Zolghadr等人17研究了CO2-烷烃体系特性随温度的变化情况。油气混相过程及界面特性会受到各种条件的影响，导致其特性有很大的差异性，但目前很少有学者从热力学角度以及分子的微观角度出发来研究不同气体造成驱油效果不同的根本原因，进而得到统一的认识。在工业驱油过程中，有一项很重要的参数—最小混相压力(MMP，Minimum Miscibility Pressure)，MMP是油气两相界面张力为0时所对应的最小气相压力26,27。当压力达到MMP时，油气两相混合均匀，界面消失28。MMP是油气达到混相的界限压力，充分认识MMP的变化规律及其微观机制对于气驱油过程的控制及优化都非常关键。

本文针对于吉林某油田的实际油组分，通过MD模拟研究CO2、N2及其1 : 1混合气与油相的混相过程，获得油气密度分布、界面张力及其MMP，揭示了不同气体组分对油气混相过程的影响规律，并从油气分子间相互作用的角度阐述了其微观机制。研究结果对深入认识气驱油混相过程中油气两相界面传质过程具有重要的意义，并可指导气驱提高原油采收率技术的优化与设计，工程价值突出。

2 MD模型

2.1 模拟系统

模拟系统中，油的组分是根据吉林某油田的实际组分模化的，通过烷烃链长来划分出短链、中长链和长链，并选取各个区域中的某一种烷烃来代替整个区域，如用C7来代替C7-C10。本模拟中的油组分具体见表1，系统温度取为地层温度(371.15 K)。

表1 模型油组分Table 1 The major component of oil in a simulation system.

构建初始模型时，在特定大小盒子里随机生成气体分子和烷烃分子，模型总尺寸为8 nm × 8 nm ×19.5 nm，采用三明治模型以保证体系的对称性12,24，沿着z方向依次为气相、油相、气相，如图1所示。本模拟研究了不同气体压力下的油气混相过程，先利用Materials Studio软件的不定型建模，在定温度371.15 K、定体积确定密度和压力的基础上确定分子的数目，然后再将盒子按照“气-油-气”叠加起来。表2展示了不同压力下的气体分子数目(仅针对一个气相盒子，整个体系的气体分子数是该分子数目的2倍)。MD模拟计算采用LAMMPS (大型原子/分子大规模并行模拟器)软件来开展。x-y方向采用周期性边界条件，z方向为反射边界条件。气体分子的数目通过371.15 K下不同压力来确定，以此模拟采油过程中向油藏中注入不同压力的气体。二氧化碳分子采用全原子EPM2模型，C－O键键长为0.116 nm，键角为180°12；N2分子采用全原子模型，N－N键键长为0.1112 nm，键角为180°28；烷烃分子采用联合原子模型(United-atom model)，把CH3和CH2等效成一个原子，利用这种模型在不影响模拟结果精度的基础上可以大大降低计算量29。键角采用Harmonic来描述其相互作用，二面角则选取opls模型来表示。采用PPPM方法来考虑长程静电力，通过带有极性项的12-6 Lennard-Jones势能模型充分考虑了范德华力和库仑力，如式(1)所示。作用力的截断半径取为10 Å (1 Å = 0.1 nm)。用Lorentz-Berthelot混合法则30来计算体系中交叉原子间的Lennard-Jones势能参数，如式(2)所示。各个原子的势能参数及所带电荷值见表3。

表2 不同压力下的气体分子数目Table 2 The number of gas molecules with different pressure.

表3 原子的Lennard-Jones势能参数以及电荷值Table 1 Lennard-Jones potential parameters of molecules and atomic charges.

图1 模拟初始构型图Fig. 1 Initial configuration of simulation system.

式(1)中，ε、σ、qi、qj、C、χ、rcut分别代表了能量参数、长度参数、原子i和原子j的电荷、静电常数、介电常数、范德华力及短程静电力的截断半径；式(2)中，ε和σ分别代表了Lennard-Jones模型中的能量参数和长度参数。

在整个模拟过程中，为了保持体系的温度体积保持不变，采用NVT系综，利用Nose-Hoover热浴维持温度，防止不同工况下油相的密度发生改变。时间步长选用1 fs，每隔10000步输出一次原子坐标信息，总模拟时长为20 ns。模拟系统在10 ns左右可达到稳定状态，形成稳定的界面。本文对最后5 ns的数据进行统计计算与分析，得到油相与气相的密度分布、界面位置以及界面厚度，然后计算不同工况下的界面张力，进而确定不同气体对应的MMP。

3 结果与讨论

3.1 密度分布

模拟每隔5 ns进行密度分布统计，发现体系在10 ns时开始达到稳定，15-20 ns密度分布曲线波动不大，并与10-15 ns统计出来的密度分布曲线对比发现几乎没有差异。因此取15-20 ns数据，计算分析出了气体分子和油分子沿垂直于界面方向(z方向)的密度分布。图2a，b分析了不同压力的纯CO2驱油下油气两相密度(ρ)在z方向(Z)上的变化情况。不难看出对于CO2气体驱替的情况，有一部分CO2分子会在油气界面处形成一个吸附层，吸附层处CO2密度达到最高。随着气体压力的升高，CO2密度增大，而油的主相密度随之减小，说明在驱油过程中油分子由于驱替气体压力的升高而膨胀，油气两相混合更均匀，油分子被更好的驱替出来。这一变化与Zhang等人22得到的CO2气体与辛烷相互扩散规律是吻合的。同时，在较低压力情况下系统达到稳定时，仅有一小部分油分子与气体混合，大部分油与CO2仍有非常明显的界面存在；当压力较高，油气两相混合程度已经很高，界面相对比较模糊(图2c-f)，说明此时的气驱油效果非常好。

图2 纯CO2气驱油下油与CO2的密度分布曲线与云图Fig. 2 The density distribution of oil and gas phase in CO2-oil system.

图3和图4分别给出了CO2与N2以1 : 1比例混合和纯N2驱油情况下油气的密度分布曲线。可以看出，随着压力的升高油气两相的密度变化和纯CO2驱油是一致的，气体压力升高导致气相密度增大而油相密度减小，说明在驱油过程中油分子由于驱替气体压力的升高而膨胀，气液两相混合更均匀，油分子被更好的驱替出来。对于1 : 1混合气，CO2的密度曲线在油气界面处密度出现峰值，而N2没有明显的峰值出现，这是CO2在界面处受到强的吸引作用力造成的，吸引力的强弱将在下文的分子平均作用力势(PMF，Potential of mean force)分析中得到具体解释。

图3 CO2与N2 1 : 1混合气驱油下油与气相的密度分布曲线Fig. 3 The density distribution of oil and gas phase in CO2/N2-oil system.

图4 纯N2气驱油下油与气相的密度分布曲线Fig. 4 The density distribution of oil and gas phase in N2-oil system.

3.2 界面厚度

根据上文分析得到的油气两相的密度分布曲线，以及吉布斯切面的定义(油的密度为其主相密度的10%-90%)可确定出油气界面位置及界面厚度(δ)，如图5a所示；对于不同的驱替气体，可统计出界面厚度与压力的关系，见图5b。对于过高的气体压力未给出界面厚度，因为此时气相与油相几乎混合均匀，界面已经变得很模糊，无法进行界定。从图5b可以发现随着压力的升高，界面厚度随之增大，这是因为更多的气体分子从气体主相扩散到油相中，与此同时油相中向气相移动的分子也增多，油气两相混合程度增强，界面变得更宽。极端情况便是超过MMP时，油气两相完全混合均匀，界面不存在。另外可以看出，当驱替气体处于同一压力状态时，随着CO2组分的增加，界面变得更厚。这是因为同等气体压力下CO2的驱替效果好于N2，油气的混合程度更明显，界面更宽。但是由于不同气体的MMP不同，他们达到接近混相状态压力范围也不一样；例如CO2在较小的压力下就能混相，所以其界面厚度在较小压力下就会变得很大。在接近混相的状态下，不同气体的界面厚度都差不多。

图5 界面示意图以及界面厚度与气体压力的关系Fig.5 Schematic diagram of interface and relationship between interface thickness and gas pressure.

3.3 界面张力与MMP

油气界面的界面张力使用吉布斯公式可以计算得出，基于压力张量的界面张力(δ)计算P公式如下24,31：

式(3)中，PN和PT分别代表法向与切向的压力，Pii(i=x,y,z)是不同方向上的压力张量，Lz为模拟系统中盒子沿着z方向上的尺寸；式(4)中，k为原子编号，ν为该原子的速度，N为体系中所有原子的数量，m为原子的质量，r为原子不同时刻的位移，f为原子受到的力，V为整个体系的体积。

通过对最后5 ns得到的Pii取值并进行平均运算，结合计算γ的公式，得到了不同压力状态下平衡后5 ns的平均界面张力，如图6所示。在计算过程中没有在界面张力接近于零的压力区间进行取值计算，是因为这个时候油气基本达到混相状态，统计出来的参数精度不高。从图6可以看出随着压力的升高，界面张力随之降低，这是因为驱替气体压力升高，油相膨胀密度降低，气相与油相的混合程度更强，故而界面张力随之减小，这与Zolghadr等人17利用实验测得的CO2与烷烃混相的规律是一致的。除此之外，Lara等人8与Ayirala等人26的结论也符合该规律。理论上，达到MMP时油气界面的张力为零；在实际计算过程中，无法统计出完全混相时的界面张力。Orr等人32明确指出利用油气界面消失可以精确地计算MMP。因此，本文通过已经确定出的不同压力时的界面张力，进行线性拟合并延伸17,29,33，得到零张力对应的压力，即为MMP33。图中有些点相对发散，这是由MD模拟中初始时刻随机速度引起的误差造成的，但这对MMP的计算结果影响不大。通过上述方法得到不同气体对应的MMP，如图7所示。可以看出，纯CO2驱替时MMP为23.0 MPa，CO2与N2摩尔比为1 : 1时得到的MMP为50.7 MPa，纯N2驱替油时获得的MMP为119.0 MPa。吉林某油田实际测得的CO2MMP为22.3 MPa，模拟值与实际值吻合得很好，说明通过MD模拟完全可以精确地描述油气混相过程及预测MMP。此外，从图7中可以看出纯CO2下的拟合线斜率最大，这说明利用纯CO2驱油并逐渐升高气相压力时体系能够很快达到油气混相状态，而另外两条线明显比纯CO2驱油下的斜率小，意味着当体系中N2比例升高后气相需要增加更大的压力才能使油气两相混合均匀。该差异体现了不同气体种类对油气混合程度的影响，为了解释这一现象，下文将对不同分子与油相之间的作用力进行详细分析与解释。

图6 不同气体压力对界面张力以及MMP的影响Fig. 6 Effects of gas pressure on interfacial tension and MMPs.

图7 不同气体驱油下的MMPFig. 7 Effects of gas components on MMP.

3.4 油气混相分子机制

从上文中可以得知，同一气体在不同压力下驱油效果不同，不同气体驱油得到的MMP也存在一定的差异。为了探究影响气驱油过程中油气混相的分子机制，首先以N2驱油为例，计算分析了不同压力下整个系统的总能量(E)随着时间(T)变化情况，如图8a所示。可以看出，随着气体压力的升高，系统总能量逐渐升高，油气混合程度高。这是因为压力越高，更多的气体分子与油分子相互作用碰撞，系统的势能与动能升高，说明油气混合程度高。这是因为压力越高，更多的气体分子与油分子相互作用碰撞，系统的势能与动能升高，油气混合程度加剧。这是针对同一气体在不同压力下，从系统能量角度对油气混合机制进行了分析。

图8 N2驱油的能量变化曲线与不同气体的PMF分布曲线Fig. 8 The time-varying energy profiles in N2 flooding and PMF profiles of different gas molecules.

为了探究CO2和N2两种不同气体驱油效果显著差异的原因，本文从微观作用力角度对这两种气体分子与油相之间的相互作用能进行了分析。利用PMF分布曲线来研究相互作用能的差异。PMF是利用LAMMPS中通过柔性谐波抑制偏置的Adaptive Biasing Force方法来进行采样计算的。计算过程中固定油相，气体分子依然按照NVT系综进行运动；沿着模拟系统的z方向划分网格，网格间距为0.02 nm，每个网格空间中的样本数设置为40。图8b所示为系统归一化高度z/L在0.1-0.9范围内的PMF分布曲线。从图中可以清楚看出，在油相所在的位置处，CO2和N2的PMF分布曲线都表现出了吸引力，但CO2的凹陷程度要高于N2，表明了油相对CO2的吸引力明显大于N2。所以在驱油过程中，同一压力下CO2与油相的混合程度更高。要获得同等的驱油效果，N2分子需要的压力要远高于CO2，N2驱油的MMP要远大于CO2。这就从微观作用力角度更深入地解释了油气的混合机制。不难看出，油气混合程度的强弱受到了体系中总能量以及气体分子与油相作用强度的影响，如图9所示。当整个体系处于较高的能量并且驱替相与油相的吸引作用较强时，油气分子向另一相扩散得更多，此时油相膨胀程度增加，整个系统更容易混合。

图9 气驱油过程中油气混相分子机制示意图Fig. 9 Schematic diagram of molecular mechanism of miscibility process in gas flooding(green parts represent miscible areas).

4 结论

本文采用MD模拟手段对吉林某油田气驱油过程进行了模拟计算，研究了不同驱替气体种类以及压力工况下，油气两相的密度分布、界面特性以及MMP，并对气驱油过程中油气混相的分子机制进行了分析。结果表明，随着驱替气体压力的升高，油气两相的界面厚度增加，油气两相的混合程度增加。值得注意的是，当气体压力足够高时，两相的混合程度变得很高，这时没有清晰的界面存在。纯CO2驱油的MMP远远小于纯N2驱油，当这两种气体以摩尔比为1 : 1混合时，得到的MMP介于两种纯气体驱油之间，说明在相同压力下CO2组分越多油气混合越容易。最后，本文通过对系统的总能量和不同气体分子PMF的分析解释了油气混相的分子机制。同一气体在压力升高时，系统的总能量随之升高，油气分子间相互碰撞增强，油气两相混合程度增加；对于不同气体分子，油相对CO2分子的吸引强而对N2的吸引弱，导致CO2对应的MMP较小。总之，本文从分子层面上获得的气驱油油气混相过程的界面传质特性及分子机制对气驱油技术具有重要理论指导。