纳米孔隙储集层动态渗吸数学模型

田伟兵，吴克柳，陈掌星, ，雷征东，高艳玲，李靖

（1. 中国石油大学（北京）油气资源与探测国家重点实验室，北京 102249；2. 加拿大卡尔加里大学石油及化学工程系，阿尔伯塔，T2N 1N4；3. 中国石油勘探开发研究院，北京 100083）

0 引言

致密油储集层孔隙以纳米尺度为主，流体的赋存方式和流动机理均与常规油藏不同[1]。目前压裂开发与注水吞吐是致密油的主要开发方式。开发过程中，在毛细管力的作用下，裂缝中的水（压裂液）经渗吸进入基质排出其中的油[2-5]。

纳米尺度孔隙的渗吸参数主要受纳米限域效应的影响。Wu等[6]研究了孔隙中纳米限域效应对水的黏度与滑脱效应的影响，得到了纳米尺度单孔单相流动条件下的泊肃叶公式。Supple等[7]开展了纳米尺度下孔隙中液驱气渗吸实验，认为渗吸不遵循 LWR（Lucas-Washburn-Rideal）方程[8]，且渗吸距离与时间是线性关系。Feng等[9]对纳米尺度孔隙中水-气界面张力的变化进行了探讨，但没有将其应用到纳米尺度渗吸过程中。Fernández-Toledano等[10]验证了纳米尺度孔隙中分子动力学理论[11]在固-气-液渗吸体系中的适用性。Kelly等[12]研究了纳米尺度孔隙中的液驱气渗吸问题，发现由于流体黏度改变、液膜等的影响，渗吸距离不能用传统的LWR方程计算。Gruener等[13]通过研究岩心中纳米限域效应对液驱气渗吸的影响，认识到弯液面的连续性受孔隙长短的影响。Liu等[14]在研究致密岩心中水驱油渗吸时发现渗吸初期采出程度与时间呈线性关系。可以看到，目前对纳米尺度孔隙中的油-水渗吸研究较少，对纳米尺度孔隙中动态接触角对水驱油渗吸的影响机理不明确。

本文在考虑动态接触角效应、纳米限域效应（包括多层黏附效应和滑脱效应）、惯性效应和入口端效应的基础上，建立了纳米尺度孔隙中油-水渗吸方程，推导了固-油-水三相接触线摩擦系数与界面区流体黏度的关系式。同时通过毛细管束模型和对数正态分布理论，得到了致密岩心中的渗吸模型。最后通过孔隙、致密岩心的渗吸模拟对模型的适应性进行了验证，并对影响渗吸的关键参数进行了分析。

1 纳米孔隙动态渗吸数学模型

1.1 孔隙中的油-水渗吸方程

基于 LWR理论，传统的孔隙油-水渗吸方程（标记为LWR模型）[8]如下：

考虑纳米限域效应、入口端效应和惯性效应，则渗吸方程（标记为ND模型）可表示为[15]：

ND模型忽略动态接触角效应（即假设接触角不随渗吸速度改变而改变），但事实上，由于接触线处具有摩擦效应，接触角将随着渗吸速度的变化而改变，是动态变化的。如考虑动态接触角效应，则纳米孔隙渗吸方程（标记为CD模型）可表示为：

（3）式中惯性力为：

基于平衡接触角计算的毛细管力为[15]：

基于考虑滑脱的泊肃叶方程，黏性力可表示为[6]：

滑脱长度可表示为[6,16]：

入口端效应可以表示为[17]：

重力作用为[15]：

三相接触线摩擦力为[15]：

（12）式中三相接触线摩擦系数表征了动态接触角效应的强弱。

1.2 三相接触线摩擦系数

根据分子动力学理论，动态接触角和速度的关系为[18-19]：

考虑到三相接触区中固相、润湿相和非润湿相分子以混合形式存在，在此提出权重法计算润湿活化自由能：

式中权重系数为单相流体黏度占总黏度的比值。

黏度效应的活化自由能通过Eyring理论表征[20]：

岩石-流体相互作用的活化自由能为[21-22]：

由（14）—（19）式可知，三相接触线的摩擦系数为：

1.3 接触线摩擦系数与固-液界面区流体黏度的关系

考虑多层黏附效应，固-液界面区黏度为：

结合（18）、（19）式，界面区流体黏度为：

综合（20）、（22）和（23）式，可得接触线摩擦系数和界面区域流体黏度的关系：

此外，流体有效黏度可通过面积权重法计算[6]。

1.4 致密岩心油-水渗吸模型

Wang等[23]研究了纳米孔隙介质中分子相互作用、动态接触角效应与入口端效应对水渗吸动态的影响，获得了较为可靠的渗吸动态数据（经过 Lattice Boltzmann验证的模拟数据）。为讨论多层黏附、滑脱、入口端、惯性等效应对渗吸动态的影响程度，在此采用该文献中渗吸模型参数作为本文 CD模型的基础参数（见表1），其渗吸动态数据作为验证数据，模拟分析忽略多层黏附、滑脱、入口端、惯性等效应对渗吸结果的影响。

表1 孔隙渗吸模型验证所需基础参数

其中有效界面张力计算公式为[24]：

界面区域的厚度取两层流体分子厚度之和[6]，润湿相界面区域密度计算公式为[25]：

假设润湿相界面区域密度与流体常规密度的比值和非润湿相界面区域密度与流体常规密度的比值相同，非润湿相界面区域密度同样采用（26）式计算。

有效密度采用权重法计算：

CD模型忽略多层黏附、滑脱、入口端、惯性等效应时，模拟结果与验证数据的对比曲线如图 1所示。可以看到，当忽略入口端与惯性效应时，模拟曲线与验证数据曲线基本重合，说明入口端与惯性效应对渗吸的影响很小，可以忽略。而忽略多层黏附效应时，模拟曲线偏离验证数据曲线较远，误差较大，该项不能忽略。

图1 CD模型忽略多层黏附等效应时的模拟结果

经计算，亲水性储集层的水滑脱长度较小，为0.10～0.41 nm；油滑脱长度与孔隙壁面性质、油的组分相关，计算参数中油的主要组分为正辛烷，油滑脱长度亦较小，为0.19 nm，该种情况下滑脱效应对渗吸动态的影响可以忽略。但是当油组分改变时，油滑脱长度将发生改变，如已知表面能为61.2 mN/m的砂岩中，正十二烷烃和正十四烷烃在其表面的滑脱长度分别为5.7 nm和28.0 nm[16]，此时滑脱效应对渗吸的影响将较为显著，该种情况下滑脱效应对渗吸动态的影响不能忽略。

基于上述分析，如忽略入口端与惯性效应的影响，此时由（3）、（5）—（7）和（12）式整理得：

为方便表征和计算，上式可表示为：

则渗吸距离和体积分别为：

考虑初始含水饱和度和渗吸后残余油饱和度的影响，岩心中可发生流体流动的孔隙体积为：

而基于毛细管束的岩心可流动孔隙体积为：

结合（32）和（33）式，岩心中可发生流体流动的毛细管总数量为：

则岩心渗吸体积、油相采出程度可表示为：

岩心孔隙尺寸可通过对数正态分布模型表征[26]：

2 模型可靠性验证

2.1 模型验证计算流程

采用表 1中的基础数据，验证动态接触角模型、孔隙渗吸模型的可靠性，具体计算流程如图2所示。

图2 模型验证计算流程

2.2 纳米尺度动态接触角模型验证

Stroberg等[27]首次开展了纳米尺度渗吸分子动态模拟研究，选取其研究结果对本文模型进行验证。首先计算界面区域的流体黏度，然后计算接触线摩擦系数（其中分子跳跃长度作为拟合参数），最后通过（13）式计算动态接触角与速度的关系。由图 3可见，本文提出的动态接触角模型与文献模型（分子动态模型）模拟结果吻合度较高（模拟孔隙半径分别为2.5 nm、4.0 nm），说明本文模型可以较好地表征渗吸过程中接触角的动态变化。

图3 纳米尺度动态接触角模型验证

2.3 孔隙渗吸模型验证

采用表 1中的基础数据，分别采用 LWR、ND、CD模型模拟渗吸时间与渗吸距离的关系，并将模拟结果与验证数据曲线进行对比（见图 4）。可以看到，传统LWR模型未考虑纳米限域效应、入口端效应、惯性效应、动态接触角效应，模拟曲线离验证数据曲线最远，误差最大，不能表征纳米尺度孔隙的渗吸过程；ND模型虽然考虑了纳米限域、入口端、惯性等效应，但未考虑动态接触角效应，其模拟曲线虽然比LWR模型更接近验证数据曲线，但误差仍然很大；在 ND模型基础上，CD模型修正了动态接触角效应的影响，模拟效果较好，拟合优度为 0.99。可见，动态接触角效应对渗吸的影响最为显著，纳米限域效应（多层黏附效应和滑脱效应）次之，惯性效应和入口端效应影响最小。

图4 孔隙渗吸模型模拟结果对比曲线

2.4 模拟渗吸采出程度验证

Ren等[28]开展了渗透率和黏度对致密岩心渗吸的影响研究，结果表明采出程度随渗透率增大而升高，随油相黏度增大而减小。王秀宇等[29]研究了温度、压力对致密岩心渗吸的影响，证实温度越高，采出程度越高，而压力存在一个最佳范围。Kathel等[30]探讨了润湿性对致密岩心渗吸的影响，发现表面活性剂可以使岩心发生润湿反转，自发渗吸采出程度约为 68%。上述文献共选用 5块实验岩心，岩心编号分别为岩心A-5[28]、岩心 R4-2[29]、岩心 R4-4[29]、岩心 2[30]、岩心7[30]，均为致密砂岩，岩心基础参数如表2所示[28-30]，通过实验获取了岩心渗吸采出程度随渗吸时间变化的5组数据。

表2 文献中实验岩心基础参数

根据表 2中的基础数据，采用 CD模型模拟渗吸采出程度与渗吸时间的关系曲线，并与 5组实验数据进行对比（见图 5）。由图可知，CD模型模拟结果与文献实验结果吻合较好，可较好地模拟致密岩心中的油-水渗吸过程。但CD模型忽略了重力影响，因此当毛细管力较弱时，重力作用对渗吸动态存在一定影响，导致模拟曲线与实验数据存在一定的偏差。

图5 渗吸采出程度与渗吸时间的关系曲线

3 渗吸动态主要影响参数

根据表1中的基础数据，采用ND、CD模型重点模拟孔隙半径、界面张力的变化对渗吸动态的影响。

3.1 孔隙半径

渗吸过程中的主要阻力为惯性力、接触线摩擦力、湿相（水）黏滞力和非湿相（油）黏滞力，其中单项阻力与这 4项阻力之和的比定义为单项阻力比。不同半径孔隙渗吸过程中单项阻力比随渗吸时间的变化如图6所示。由图可知：①渗吸初期，以惯性力为主导，以半径为 50.0 nm的孔隙模拟为例，当渗吸时间小于0.100 ns时，随着时间增加，惯性阻力比逐渐减小，接触线摩擦阻力比逐渐增加。当渗吸时间为0.1～100.0 ns，接触线摩擦阻力比最高、油相黏滞阻力比次之，接触线摩擦力占主导地位，且接触线摩擦阻力、油相黏滞力基本保持不变；渗吸时间大于100 ns时，接触线摩擦阻力比和油相黏滞阻力比逐渐减小，而水相黏滞阻力比逐渐增加。②接触线摩擦力不仅作用于渗吸后期，呈下降趋势，也作用渗吸初期，且阻力比呈上升趋势。③随着孔隙半径减小，惯性阻力比曲线左移，惯性力作用时间缩短；渗吸初期接触线摩擦阻力比曲线向左移动（对应孔隙半径从50.0 nm到20.0 nm再到10.0 nm的模拟曲线变化。由于本文最小模拟时间为0.001 ns，因此孔隙半径为5.0 nm和2.5 nm的模拟曲线未能完整显示该过程），动态接触角效应增强，渗吸中后期接触线摩擦阻力比曲线向下移动，动态接触角效应减弱；油相黏滞阻力比曲线向上移动，油相黏滞力的影响增加；水相黏滞阻力比曲线向左移动，水相黏滞力初始作用时间提前。

图 7为不同渗吸模型在不同孔隙半径条件下渗吸时间与渗吸距离、渗吸距离增强因子（CD模型渗吸距离与 ND模型渗吸距离的比值）的关系，可以看到，孔隙半径减小，CD、ND模型的模拟曲线位置下降，相同渗吸时间内渗吸距离减小。随着渗吸时间的延长，渗吸距离增强因子曲线呈“U”形，早期下降，中期趋于平稳，后期上升，因渗吸距离与动态接触角效应负相关，故动态接触角效应表现为早期上升，中期趋于平稳，后期下降，与图6b中的接触线摩擦阻力比变化一致。同时还可以看到，孔隙半径减小，渗吸初期曲线向左移动，动态接触角效应的影响增强，渗吸中后期曲线向上移动，动态接触角效应的影响减弱，仍与图6b中的接触线摩擦阻力比变化一致。

图6 不同孔隙半径条件下渗吸时间与单项阻力比的关系

图7 不同渗吸模型在不同孔隙半径条件下渗吸时间与渗吸距离、渗吸距离增强因子关系曲线

图 8为不同孔隙半径下，渗吸时间与动态接触角的关系曲线。渗吸初始，水相刚接触孔隙内壁面时渗吸速度近似为零，根据（13）式可知，此时的接触角为平衡接触角（初始时间接近于零，远小于本文最小模拟时间0.001 ns，因此图中未显示该过程）；随着渗吸时间增加，渗吸速度增大，接触角逐渐增加至最大，并在一定时间内保持稳定，如图 8所示；随着渗吸时间进一步增加，渗吸速度减小，动态接触角逐渐减小，如图 8所示。根据（13）式可以预计，当渗吸时间增大至渗吸速度减小趋近于零时，接触角将最终趋近平衡接触角（趋近时间远大于本文最大模拟时间1 000 ns）。

图8 不同孔隙半径下渗吸时间与动态接触角的关系曲线

3.2 界面张力

油水界面张力越大，渗吸动力越大，然而流体与岩石之间的黏附功也会增加。图 9为不同界面张力条件下渗吸时间与接触线摩擦阻力比的关系曲线，可以看到，界面张力增加，曲线上移，接触线摩擦阻力比增加，动态接触角效应增强。

图9 不同界面张力下渗吸时间与接触线摩擦阻力比的关系曲线

图 10为不同界面张力条件下渗吸时间与渗吸距离、渗吸距离增强因子的关系曲线。可以看到，界面张力由0.01 N/m增大至0.05 N/m，CD、ND模型的模拟曲线位置上升，相同渗吸时间内渗吸距离增加；但当界面张力由0.05 N/m增大至0.10 N/m时，CD模型的模拟曲线位置下降，ND模型的模拟曲线先上升后下降，可见界面张力存在临界值（本算例位于 0.05～0.10 N/m）。因此，致密储集层渗吸提高采收率，过低的界面张力并不能获得更好的渗吸效果，而应该通过优化获得最佳的界面张力。图 10c渗吸后期的结果同样显示界面张力存在临界值。

图10 不同界面张力条件下渗吸时间与渗吸距离、渗吸距离增强因子的关系曲线

图11为采用CD模型模拟的不同界面张力条件下渗吸时间与动态接触角的关系曲线。随着渗吸时间的延长，动态接触角从最初的平衡接触角逐渐增加至最大，并在一定时间内保持稳定，随后动态接触角将逐渐减小并趋近于平衡接触角。同时最大动态接触角随着界面张力增加而增大，即界面张力增加，接触角的变化范围增大，动态接触角效应增强。

图11 界面张力对动态接触角的影响

4 结论

纳米孔隙渗吸过程中，动态接触角效应对渗吸的影响最为显著，纳米限域效应（多层黏附效应和滑脱效应）次之，惯性效应和入口端效应影响最小。

渗吸初期，惯性力的作用减小，接触线摩擦力的作用增加，动态接触角从初始的平衡接触角逐渐增大至最大并基本保持稳定；渗吸后期，接触线摩擦力的作用减小，接触角则从最大动态接触角逐渐降低并趋近于初始的平衡接触角。

孔隙半径减小，渗吸初期动态接触角效应增强，渗吸中后期动态接触角效应减弱；油水界面张力增加，渗吸动力增大，动态接触角效应增强。

界面张力对渗吸动态的影响存在临界值，致密储集层渗吸提高采收率，过低的界面张力并不能获得更好的渗吸效果，最佳的界面张力需通过优化获得。

本文模型考虑了动态接触角效应、纳米限域效应、惯性效应和入口端效应的影响，可较好地描述纳米孔隙和致密储集层中的渗吸过程。

符号注释：

A1——中间变量，Pa·s/m；Ai——界面区横截面积，m2；At——总流动区横截面积，m2；B——中间变量，Pa·s；C——中间变量，N/m；Cf1——与孔隙半径对数标准差相关的常数，m；Csw——与水相滑脱效应相关的常数，m；d——孔隙直径，m；Dc——岩心直径，m；f1，f2——孔隙半径的对数标准差和对数均值，无因次；f（R）——孔隙半径分布频率，m-1；f(θeq)——与固-液界面区流体密度和平衡接触角相关的函数，kg/m3；g——重力加速度，m/s2；h——普朗克常数，6.626×10-34J·s；kB——玻尔兹曼常数，1.38×10-23J/K；K——跳跃频率，Hz；ls——滑脱长度，m；L——孔隙长度，m；Lc——岩心长度，m；n——单位面积上的吸附点位个数，m-2；NA——阿伏伽德罗常数，6.02×1023；Nm——岩心中可发生流体流动的毛细管数量；p——压力，Pa；Q——单根毛细管渗吸体积，m3；Qc——岩心渗吸体积，m3；R——孔隙半径，m；Rave——平均孔隙半径，m；Ro——采出程度，%；S——相邻分子可进入孔隙的等效壁面积，m2；Swc——束缚水饱和度，%；Sor——残余油饱和度，%；t——渗吸时间，s；T——温度，K；v——速度，m/s；V——单个流体分子的有效体积，m3；Vm——岩心中可发生流体流动的孔隙体积，m3；x——渗吸距离，m；x（R,t）——渗吸距离函数，m；ΔG——活化自由能，J；α——由固体构成的等效壁面积比例，无因次；γ——界面张力，N/m；δ——相邻层油分子间的平均距离，m；δl——Tolman长度，m；ζ——接触线摩擦系数，Pa·s；θd——动态接触角，（°）；θeq——平衡接触角，（°）；λ——跳跃长度，m；μ——黏度，Pa·s；ρ——密度，kg/m3；σ1——真空中油相的表面张力，N/m；τ——迂曲度，无因次；τt——弛豫时间，s；φ——孔隙倾角，（°）；φ——孔隙度，%；ω——权重系数，无因次。下标：b——常规条件下流体参数值；ce——基于平衡接触角计算；e——流体参数有效值；en——入口端；fl——接触线摩擦；g——重力；i——界面区；in——惯性；l——液相；max——最大值；min——最小值；nw——非润湿相；s——固相；v——黏性；w——润湿相；wet——润湿。