页岩油微尺度流动表征及模拟研究进展

苏玉亮，王 瀚，詹世远, 2，王文东，徐纪龙

1)中国石油大学(华东)石油工程学院， 山东青岛 266580；2)成都理工大学能源学院，四川成都 610059

随着水平井和多级压裂技术的发展，页岩油效益开采成为一种非常重要的非常规石油资源[1]．美国作为世界页岩油产量最高的国家，2019年页岩油产量增加到28.3亿桶(标准桶，standard tank barrel, STB)，占全球总石油资源产量的63%[2]．2020年中国石油进口依存度高达74%，远超过50%的安全警戒线，能源安全供应风险亟需高度关注，加快页岩油勘探开发是保障国家石油能源安全的途径之一[3]．

北美页岩油资源有效开采率先取得了重大突破，引发了世界范围内的能源革命[1]．北美页岩油产层主要分布于海相或前陆盆地，具有沉积面积大、连续性好、压力系数高和成熟度高等特点．中国页岩油主要产层多分布于陆相沉积的坳陷和断陷型盆地，具有非均质性强、压力系数整体偏小和成熟度低等特点[4]．页岩储层普遍具有超低渗透率(1×10-9～1×10-4μm2)，孔隙直径为纳米级(1～200 nm)[5]，孔隙结构复杂[6]．在受限纳米孔隙中，固-液分子间作用力造成复杂的流体性质和壁面液固流动性质，此时传统的宏观达西渗流方程不再适用于准确表征纳米受限空间中的流体流动特性[7-8]．目前，针对页岩油-气-水多相流体在纳微空间中复杂流动行为，学者已经进行了大量研究[9-13]．

本研究旨在总结和评述近年来页岩油微观流动机理理论、分子动力学模拟和孔隙尺度流动模拟方面的进展和未来研究趋势．首先，介绍页岩油在受限空间中受固-液分子间作用力影响的复杂流体性质和壁面流动性质；并基于这些性质总结4种单纳米孔隙中流体流动表征方法；最后介绍基于分子动力学模拟(molecular dynamics simulation, MDS)、格子Boltzmann方法(lattice Boltzmann simulation method,LBM)、孔隙网络模型(pore network model, PNM)和LBM-PNM耦合方法等模拟页岩油流动过程，明晰页岩油流动机制研究现状．

1 页岩油微观渗流机理及理论表征

1.1 壁面滑移和流体性质非均质性

受液-固之间相互作用力的影响，流体在固体表面真实速度一般大于0，该速度被定义为滑移速度．滑移速度可通过壁面剪切速率和滑移长度计算得到．图1(a)为纳米孔隙横截面，图1(b)为不同滑移速度和黏度非均质性下的速度剖面，滑移速度vs由0逐渐增大;ls为滑移长度，单位：nm．速度剖面中实线为基于有效黏度计算的速度剖面，虚线为基于近壁面相油黏度计算的速度剖面，该速度剖面用于计算壁面滑移速度．纳米孔隙中页岩油被分为近壁面油和体相油相，近壁面油的厚度为δ， 一般为2个分子层厚度，油相厚度约为0.98 nm[14]，水相厚度约为0.7 nm[15-16]．多孔介质中孔隙结构复杂多样，拐角处存在壁面作用力叠加影响，此时近壁面流体的区域大小将会发生变化，可以通过计算准确的密度分布来进一步确定多孔介质中的近壁面区域．

图1 不同水相接触角下速度剖面示意图[8]Fig.1 (Color online) Velocity profiles with different solid-liquid molecular interaction forces[8]

对于油相来说，WANG等[14]采用MDS方法研究了油在直径为5.24 nm的石英和石墨孔隙中的滑移流动，模拟条件下油在亲水石英孔隙中的滑移长度远小于亲油石墨孔隙中滑移长度，该研究结果与水相理论公式和常规认识(润湿性越强，滑移长度越小)相悖，原因是分子尺度上石英表面较为粗糙，而石墨孔隙表面光滑，这也间接说明滑移长度受表面粗糙度影响较大．WANG等[9, 20-21]研究了页岩油在亲油有机质孔隙和亲水无机质孔隙中的流动，但由于滑移参数的选择并不能真实反应实际地层纳米孔隙中的状态，只定性分析了不同滑移参数、孔隙几何形状条件下页岩的油流动规律．流体在孔隙壁面的滑移流动机制受许多物理参数的影响，包括壁面非均质润湿性[19, 22]、壁面粗糙度[23-24]、剪切速率[25]、纳米气泡或气膜[26]、流体极性[27-28]、流体黏度[29]、温度[30-31]、孔隙尺寸[32]、压力梯度[33]和烃链长度[30]等．FALK等[34]在直径分别为1 nm和7 nm的孔隙中，发现滑移长度可高达500 nm和120 nm，这比MAJUMDER等[35]的测量结果小很多．一般来说，MDS中采用的孔隙具有理想的光滑壁面，然而，在实验过程中多因素影响下测量的滑移长度约为MDS计算滑移长度的14.6倍[7]．页岩储层中纳米孔隙、流体性质和地层温压条件复杂多样，在多物理因素影响下的滑移参数复杂多样，难以定量表征．滑移长度的影响一般通过增强因子进行表征，流动增强因子定义为纳米孔隙中实际流体体积通量与无滑移Hagen-Poiseuille(HP)方程计算的体积通量之比．WU等[7-8]的理论研究表明，在壁面润湿性不同的情况下，与经典的无滑移HP方程预测的体积流量相比，液体在纳米孔隙中的流动能力提高了1～1×107倍，即增强因子可达到1～1×107．

液固分子间相互作用力除了会导致壁面滑移速度外，还会导致流体密度和黏度的非均质性．在纳米孔隙中，近壁面流体的密度和黏度不同于体相流体的密度和黏度．根据MDS，近壁面区域约有2层分子受壁力分子作用力的影响，油相厚度约为1 nm[36]，水相厚度约为0.7 nm[37]．对于水相来说，近壁面水的黏度可以通过接触角定量表征．WU等[7]通过拟合MDS和实验数据，得到近壁面水相黏度、体相水黏度和接触角之间的线性关系．结果表明，近壁面水相黏度随水相接触角的增大而增大．ZHANG等[38]基于MDS得到了近壁面水密度、体相水密度和接触角之间的线性关系．对于油相，WU等[30]研究了纳米限域效应对正构烷烃流动行为的影响，结果表明，正构烷烃流动特性受第1层烃分子界面阻力和其余烃分子黏滞阻力控制，并且依赖于壁面能和烃链长度的影响．在无机纳米孔中存在油水体系时，界面水-油相的黏度可以比体相水-油相的黏度小或大几倍．在不同近壁面流体黏度和密度的影响下，流动增强因子约为0.1～10．

1.2 理论表征方法

在纳米孔隙中，固-液分子间作用力导致的滑移边界和流体黏度/密度非均质性对流体流动能力的影响不可忽略，所以常规的HP流动方程不足以表征受限流体在纳米孔隙中的流动规律，需要对其进行纳米尺度效应修正．针对单纳米孔隙，目前修正模型主要有分区模型[20]、有效黏度模型[7]、表观黏度模型[8]和表观滑移长度．

流体在圆柱形孔隙中的HP流动方程可表示为

(1)

对速度方程进行积分，可得不考虑纳米尺度效应时流体通过孔隙的体积流量为

(2)

1.2.1 分区模型

由于流体黏度的分均质性，近壁面流体和体相流体可以分别使用方程进行描述．考虑滑移边界，近壁面流体和体相流体的流动方程[20]分别为

(3)

(4)

其中，vwall和vbulk分别为近壁面流体和体相流体流动速度，单位：nm/s；μwall和μbulk为壁面相流体和体相流体的黏度，单位：mPa·s；C1和C2分别为和滑移边界有关的参数；δ为近壁面流体厚度，单位：nm．

滑移边界条件和近壁面流体和体相流体的速度连续性条件为

(5)

vwall|r=R-δ=vbulk|r=R-δ

(6)

根据边界条件，式(3)和式(4)可改写为

(7)

(8)

对速度方程进行积分求和，分区模型得到考虑纳米尺度效应的体积流量为

(9)

(10)

1.2.2 有效黏度模型

分区模型采用两套方程分别表征近壁面流体和体相流体流动，较为复杂．THOMAS等[17]通过面积加权将流体性质非均质性用一个有效黏度进行表示．WU等[7]采用有效黏度和滑移速度对HP方程进行修正，

(11)

通过滑移长度和壁面剪切速率可得滑移速度为

(12)

根据式(11)和式(12)，基于有效黏度的速度方程为

(13)

通过对速度积分和有效黏度模型可计算得到体积流量为

(14)

根据达西定律，可得基于有效黏度的表观渗透率为

(15)

1.2.3 表观黏度模型

基于有效黏度模型，通过表观黏度综合考虑黏度非均质性和边界滑移，可求得r处的表观黏度速度方程[8]为

(16)

其中，μapp为表观黏度，单位：mPa·s；

通过速度积分可得表观黏度模型计算的体积流量为

(17)

上述方程中表观黏度为未知值，可通过有效黏度模型计算的体积流量来表征．设式(14)的体积流量与式(17)体积流量相等，

(18)

根据式(18)，表观黏度为

(19)

基于式(17)体积流量和达西定律，可得基于表观黏度方法的表观渗透率为

(20)

1.2.4 表观滑移长度模型

与表观黏度方法相同，采用表观滑移长度综合考虑黏度非均质性和边界滑移，可得速度方程为

(21)

其中，lsapp为表观滑移长度， 单位： nm．

通过速度积分可得体积流量为

(22)

根据体积流量相等的原则，表观滑移长度为

(23)

通过达西定律，可得基于表观滑移长度的表观渗透率为

(24)

根据分区、有效黏度、表观黏度和表观滑移长度模型的速度方程，当接触角分别为60°和120°时，速度剖面如图2[7-8,20]．从图2可见，分区模型中，速度剖面由近壁面和体相流体两个速度分布组成．在有效黏度、表观黏度和表观滑移长度模型中，通过面积加权平均将由于分子间相互作用力导致的非均质黏度考虑成固定黏度，此时单一黏度下速度剖面为典型的抛物线型．分区、有效黏度和表观滑移长度模型中存在非零边界速度，在真实流动情况下，滑移速度始终大于0．当近壁面流体黏度较体相流体黏度增大时，会出现表观滑移长度小于0的情况．表观黏度模型可将复杂的滑移流动简化成采用经典无滑移HP方程．根据达西定律计算得到表观渗透率，对4种模型的速度分布进行积分，不同模型的表观渗透率与接触角的关系如图3[7-8,20]．由图3可见，4种模型均能准确表征液体在纳米尺度空间中的流动能力．

图2 分区模型、有效黏度模型、表观黏度模型和表观滑移长度模型的速度剖面[7-8,20]Fig.2 (Color online) Velocity profiles of region-separation, effective viscosity, apparent viscosity and apparent slip length models[7-8,20]

图3 分区模型、有效黏度模型、表观黏度模型和表观滑移长度模型计算得到的表观渗透率[7-8,20]Fig.3 (Color online) Apparent permeability of region-separation, effective viscosity, apparent viscosity and apparent slip length models[7-8,20]

2 页岩油微观渗流模拟方法

虽然页岩油多孔介质中的流体流动可以通过岩心流动实验来研究，给出渗透率等参数，但实验不可避免地存在着昂贵和耗时等缺点．因此，数值模拟越来越受到关注[39]．多孔介质中流体流动模拟有两种方法，即直接模拟和孔隙网络模型[40]．其中，直接模拟法包括MDS(图4[36])和LBM(图5[41])．

图4 无机质孔隙油相流动模拟[36]Fig.4 (Color online) Simulation of oil flow in inorganic pores[36]

图5 多孔介质中不同边界滑移长度下水相无因次速度分布[41]Fig.5 (Color online) Dimensionless velocity distribution in porous media with different slip length[41]

2.1 分子动力学模拟

MDS方法目前已较为成熟，其模拟过程人为设定因素少，模拟结果可信度高，成为揭示纳米尺度流体运移机制的重要工具，并且随着计算机技术的高速发展，该方法逐渐在油气渗流领域得到了广泛应用．由于页岩储层中矿物成分多样，石英和黏土等无机质矿物与干酪根等有机质共存，多位学者对不同壁面类型的纳米孔隙内流体赋存状态、扩散能力与流动特征开展了大量研究[42-45]．对于无机质纳米孔隙，文献[42, 44]分析了烷烃流体在蒙脱石、石英及石墨狭缝孔隙中的分布特征与自扩散系数，研究表明，壁面附近区域烷烃流体密度出现震荡，呈现明显的分层现象；随着温度增加、压力降低、孔隙增大、碳链数目降低，烷烃流体的自扩散能力显著增大，同时发现扩散能力不仅与壁面矿物类型有关，而且呈现各向异性，如石英纳米孔隙内辛烷在平行于壁面的方向的扩散能力大于其在垂直于壁面方向的扩散能力．在外力驱动条件下，烷烃流体在光滑的石英无机矿物表面存在明显的速度滑移[36]，但孔隙中间的流体速度仍为抛物线形状．而在有机质孔隙研究方面，部分学者利用简单碳纳米结构，如多层石墨烯[46-48]、碳纳米管[34, 49-50]及碳纳米孔[51]等来表征有机质孔隙，发现外力驱动时烷烃流体的速度剖面趋于活塞状，其滑移长度可达90～185 nm[14, 49]．不过，也有学者指出，如此异常高的流量是由光滑平整的壁面表面决定的[52]，所采用的简单碳纳米结构难以代表储层内真实的干酪根孔隙表面[53]．基于干酪根的元素组成、化学基团构成和孔隙分布等特征[54]，有学者建立了较为逼真的三维无定型干酪根分子模型[55-56]，并用以开展流体运移机制研究．由于干酪根表面存在分子尺度的粗糙度，壁面附近并未出现明显的密度震荡区域，且在外力驱动时，未发现壁面滑移特征[57-59]，单组分的烷烃流体仍呈现抛物线特征，体相流体以黏滞阻力为主，含有沥青质等极性分子的多组分流体速度剖面则出现了速度平台；而在干酪根基质内的亚纳米孔隙中，烷烃流体为非达西流动，以壁面摩擦力为主，难以用常规的连续性介质假设进行描述，因此，OBILIGER等[60-61]建议利用自由体积理论对其进行描述．

2.2 格子Boltzmann方法和孔隙网络模型

在常规连续性方程无法给出准确纳米空间中流动行为时，MDS成为模拟流体在纳米空间最流行可靠的工具．然而，当前的MDS由于计算资源、计算难度的限制，在模拟具有复杂孔隙结构的多孔介质内的流体流动时面临较大挑战．

因LBM易于处理复杂结构下流体流动，而成为一种模拟多孔介质中流体流动的有效工具．通过并行计算和算法简化，很大程度上减少了计算资源[62-63]．LBM方法的难点主要是如何准确考虑黏度非均质和边界滑移影响．WANG等[41]基于LBM研究了水相在纳米多孔介质中的流动规律，定性讨论了润湿性影响下滑移长度和黏度非均性对流动行为的影响规律．ZHAO等[64]基于分子模拟得到的滑移长度和有效黏度，采用多松弛格子Boltzmann方法研究了水和烷烃在石英纳米多孔介质中的流动规律，并探讨了接触角和出入口效应的影响机制．ZHANG等[38]将分子间相互作用力耦合到LBM计算中，通过LBM模拟得到流体在纳米孔隙中的密度分布和速度分布．ZHAO等[65]采用格子Boltzmann方法研究了单纳米孔隙几何结构对限域水流动的影响，并基于模拟结果建立了不同形状纳米孔隙流动的经验公式．CHENG等[66]采用MRT-LBM研究了水在粗糙纳米孔隙中的流动规律，结果表明，随着表面相对粗糙度的增加，疏水粗糙纳米孔隙的输运能力显著增强，而亲水性纳米孔隙的输运能力变化不大．孔隙网络模型从真实孔隙结构中提取，基于简单的圆管模型来模拟流体在简化多孔介质中的流动规律．与直接模拟方法计算页岩油渗透率相比，PNM虽然计算准确性较低，但计算消耗资源小，计算效率快．YANG等[67]基于修正的页岩油流动方程建立了新的孔隙网络模型，分析了滑移长度、吸附相与体相黏度比和有机质含量对页岩油渗透率的影响．大多数研究中PNM结构是典型的球棍模型，忽略了真实孔隙结构的特点．随后部分学者将直接模拟方法计算准确和PNM计算效率高的优点相结合，建立了PNM和直接模拟相结合的方法．MIAO等[68]通过神经网络预测、训练得到不同横截面孔隙的无量纲水力阻力，将其应用到PNM中，可将误差控制在20%以内．ZHAO等[69]采用LBM计算流体通过真实孔隙结构的水动力学阻力，并将其应用到PNM中，该耦合方法达到了需求精度，且计算成本更低．

目前，大多数学者采用LBM研究水在纳米多孔介质中的流动规律，所建模型都可用来研究页岩油在纳米多孔介质中的流动．但是需要通过分子模拟得到页岩油在干酪根、石英和蒙脱石等纳米孔隙中的速度分布、密度分布等微观参数，并通过LBM单纳米孔隙模拟拟合速度分布和密度分布，得到滑移长度和固液相互作用力等格子参数，最后基于该拟合格子参数，将单纳米孔隙页岩油流动扩展到复杂多孔介质结构流动，进一步研究非均质润湿、孔隙几何结构等复杂参数对表观渗透率的影响机制．因此，目前最主要的问题是如何准确模拟得到页岩油在真实多孔结构中的滑移机制．

3 展 望

页岩储层既有亲油的有机质干酪根孔隙，也存在亲水的无机质孔隙．无机质孔隙中常存在一些原生水，在水力压裂开采影响下，压裂液通过自发渗吸或强制渗吸进入页岩孔隙[70-71]．随着开采压力的降低或在近开采井处，低组分烃类会析出形成气体．因此，页岩地层中多是油水/油气/油气水多相多组分流动，而相对渗透率是了解复杂多相流动的可行性工具．目前，部分学者对基于页岩单纳米孔隙两相分层流动展开了研究．ZHANG等[37]通过考虑滑移边界条件修正HP方程，提出了纳米多孔介质中气-水两相流动的相对渗透率模型．WANG等[10]在文献[37]模型的基础上，建立了考虑滑移边界和黏度非均质性的纳米孔隙介质中油水流动相对渗透率模型，将单纳米孔隙中油水两相分为近壁面水相、体相水、油水过渡相和体相油．ZHAN等[72]利用MDS研究了油水两相在纳米孔中的流动，并在文献[10]模型的基础上考虑了液-液(油水界面)滑移．基于孔隙尺度模拟方法，WANG等[70]采用LBM和数学模型研究了饱和油纳米多孔介质中自发渗吸行为．然而，上述基于单孔隙流体流动的理论研究难以满足工程应用的实际需要，尤其是考虑到两相分布的复杂性．最近，ZHANG等[73]从孔隙尺度模拟角度研究了页岩气-水两相流动，但所提模型不适应滑移边界条件．

综上所述，关于页岩油微尺度流动的研究目前仅限于定性研究．如何准确获得页岩油/气/水在页岩有机质孔隙和无机质孔隙中的滑移参数和非均质黏度/密度参数，并考虑壁面粗糙度、迂曲度和非均质润湿等因素影响；建立考虑滑移边界和黏度/密度非均质性的多相多组分流动孔隙尺度模拟方法(LBM和PNM)是目前定量研究亟需解决的问题．

结 语

本文综述了近年来纳微受限空间页岩油流动机理研究现状，包括表征单纳米孔隙流体流动的4种理论模型(分区模型、有效黏度模型、表观黏度模型和表观滑移长度模型)和模拟方法(分子动力学模拟、格子Boltzmann模拟和孔隙网络模型)．基于理论模型可知，在理想滑移边界条件下，受边界滑移和黏度/密度非均质影响，页岩油流动能力是常规连续性方程预测流动能力的0.1～1×107倍．基于单纳米孔隙理论模型，采用间接模拟方法PNM研究流体在多孔介质中的流动，具有计算效率高和耗费计算资源小的优点，但是PNM将复杂多孔介质结构简化成简单的球棍结构，忽略了复杂孔隙结构的影响，计算渗透率等参数准确性较低．直接模拟方法MDS由于计算资源和计算难度的限制，对流体在具有复杂孔隙结构的多孔介质内的流动模拟具有较大挑战性，但是可以模拟得到页岩油流动的微观参数：边界滑移和黏度/密度非均质性应用到孔隙尺度模拟中．LBM因易于处理复杂结构和易并行减少计算资源而成为一种模拟多孔介质中流体流动的有效工具，并基于数字岩心技术可以准确计算页岩油流动渗透率．

总的来说，页岩油微观流动机制研究是基于MDS得到页岩油微观运移参数(滑移长度、体相和近壁面相流体黏度和密度等)，应用到孔隙尺度模拟PNM和LBM中，但是目前还仅限于定性研究，如何将孔隙尺度模拟结果升尺度到现场应用仍是目前研究的难点．