裂缝性页岩油藏多孔介质渗吸数学模型及应用

李冠群，苏玉亮，董奇，王飞，孙庆豪，王文东

(1.中国石油大学(华东) 非常规油气开发教育部重点实验室，山东 青岛，266580；2.中国石油大学(华东) 石油工程学院，山东 青岛，266580；3.长庆油田分公司油气工艺研究院，陕西 西安，710018；4.低渗透油气田勘探开发国家工程实验室，陕西 西安，710018)

中国非常规油藏储量丰富，是未来国内原油稳产增产的重要领域[1]。在大规模水力压裂作业后的页岩油藏中，自发渗吸是提高采收率(EOR)的最重要机制之一[2]。页岩储层孔隙结构复杂，孔径极小[3]，并且天然裂缝较发育[4]，基质孔隙和微裂缝中油水界面的毛细管力非常大，定量表征压裂液在多孔介质中的渗吸至关重要[5]。

流体向多孔介质中的渗吸是一种普遍存在的自然现象[6]，早在20世纪初，国外学者就开始了渗吸驱油的理论和实验研究。LUCAS[7]在计算中将多孔介质简化为毛管束，建立了渗吸的毛管束模型。随后，WASHBURN[8]改进了该模型，形成了沿用至今的LUCAS-WASHBURN(L-W)模型，其描述了渗吸前缘和时间的关系。在L-W 模型基础上，人们考虑了黏性力、重力[9]、惯性力[10]对渗吸的影响，对模型中的流体动力和参数进行了修正。为了简化计算过程，很多情况下忽略了非润湿相的影响，尤其是油相的重力和黏性力。但是在页岩油藏中，由于油相和水相的黏度不同，黏性力对渗吸的影响不可以忽略。同时，需要考虑重力对渗吸的影响。CAI等[11]通过考虑不同尺寸和形状的孔隙、多孔介质中自吸流线的迂曲度和初始润湿相饱和度，建立了广义自吸模型。王付勇等[12]基于毛管束模型，考虑束缚水和残余油饱和度，构建岩心尺度基质渗吸驱油数学模型，并通过渗吸实验进行了验证。WU等[13]推导了油水两相渗吸的隐式解析解，利用分形理论建立了自发渗吸和驱替对致密储层开发影响的定量评价模型。但是，微裂缝中的自发渗吸机理与页岩基质中的自发渗吸存在明显差异，他们没有考虑微裂缝中饱和油的自发渗吸。CHENG等[14]首次提出了重力作用下自吸水进入单个饱和油裂缝的隐式解析解，并提出了饱和油压裂多孔介质自吸的岩心尺度数学模型。研究认为有/无重力的渗吸速度差异随着裂缝孔径的增加而增大，并且单条饱和油裂缝存在自发渗吸临界裂缝长度。WANG 等[15]推导了润湿相自发渗入气体饱和裂缝多孔介质的半解析模型，对裂缝中有无重力影响的解析解进行了研究，并提供了重力效应评估诊断图。但是，定量表征在天然裂缝发育的页岩储层中的渗吸研究依然较少。

目前尚缺乏同时考虑基质孔隙和微裂缝中自发渗吸的定量表征方法，无法明确页岩油藏体积压裂后闷井渗吸机理。因此，构建考虑基质孔隙和微裂缝的自发渗吸隐式解析解，推导岩心尺度不同孔隙孔径、裂缝开度分布、裂缝形态的渗吸数学模型，定量评价不同影响因素对渗吸的影响，对明确页岩油藏渗吸驱油规律及提高采收率具有重要意义。

1 圆形毛细管和板状管中的渗吸模型

针对页岩储层体积压裂后采取的闷井增产技术，目前已进行了许多试验和宏观数值模拟，并取得较好效果。但是，在天然裂缝发育的多孔介质中的闷井渗吸增产机理尚不清楚，需进一步研究。

基于孔隙和微裂缝分布，将渗吸流动通道分为连通孔隙通道和微裂缝通道，分别等效成平行毛管束和板状管束，建立物理模型，如图1所示。

图1 岩心中毛管束和板状管束的渗吸示意图Fig.1 Schematic diagram of imbibition of capillary tube bundle and plate tube bundle in rock core

1.1 毛细管中的渗吸模型

首先，分析单一毛细管情况下的渗吸动力学模型是后续研究多孔介质渗吸特征的基础工作之一。建立毛细管渗吸流动力学模型，假设：

1) 毛细管为水湿孔隙；

2) 忽略入口效应和损失；

3) 惯性力和渗透压忽略不计；

4) 毛细管力是恒定的，采用静态接触角。

毛细管中的水-油分布如图2所示。图中，L为毛细管长度；v为渗吸速度，m/s；θ为接触角，(°)；r为毛细管半径，m；x为渗吸距离，m。

图2 毛细管渗吸流动模型Fig.2 Capillary imbibition flow model

根据牛顿第二定律，自发渗吸动力学方程为

式中：Fc为毛细管力，N；Fv为黏性力，N；Fg为重力，N；m为毛管中流体质量，kg；a为加速度，m/s2。

根据Young-Laplace 方程，圆形毛细管的毛细压力pc为

式中：σ为界面张力，N/m。毛细管力Fc为

根据牛顿流体内摩擦定律，毛管壁上的剪切力τv为

式中：μv为流体黏度，Pa·s；vmax为直毛细管横截面的最大渗吸速度，m/s。渗吸过程为不可压缩牛顿流体在毛细管内的缓慢流动，流体流动更趋向于稳定的Poiseuille 流动，流动阻力遵循Hagen-Poiseuille规则，最大渗吸速度vmax和平均速度vˉ符合如下关系[16]：

因此，黏性力表达式如下：

式中：Fvw为水相黏性力，N；Fvo为油相黏性力，N；τvw为水相剪切力，N/m2；τvo为油相剪切力，N/m2；μw为水相黏度，Pa·s；μo为油相黏度，Pa·s；dx/dt为渗吸速度，m/s，可表示为

毛细管中水和油总质量m为

重力的表达式如下：

式中：ρw和ρo分别为水相和油相流体密度，kg/m3；g为重力加速度，取9.8 m/s2。

将式(3)、式(6)、式(7)和式(8)代入式(1)，微分后得

在毛细管中流速较小，因此，惯性力项很小，此处可以忽略不计，式(10)简化为

对式(11)进行变形、积分后可得毛细管中自发渗吸的隐式解析解tp为

式中：A=(ρo-ρw)r2g；B=2rσcosθ-r2ρogL。当dt足够小时，误差可以忽略，式(11)差分并积分后，可得毛细管自发渗吸的数值解xpi为

第i根毛细管的渗吸量(产油量)qpi为

1.2 板状管中的渗吸模型

在本文模型中，将连通微裂缝等效为截面为长方形的弯曲板状管，分析单一板状管情况下的渗吸动力学模型，与毛细管渗吸模型共同研究多孔介质的渗吸特征。建立单一板状管的渗吸流动力学模型，假设：1) 板状管为水湿；2) 忽略入口效应和损失；3) 惯性力和渗透压忽略不计；4) 采用静态接触角。单一板状管中的油-水分布如图3所示。

图3 单一板状管中自发渗吸示意图Fig.3 Schematic diagram of spontaneous imbibition in a single plate tube

根据Young-Laplace 方程，裂缝中油水界面的毛细力Fc为

式中：w为板状管开度，m；l为宽度，m。

根据式(6)计算，黏性力Fv为

重力的表达式如下：

将式(15)、式(16)和式(17)代入式(1)，微分后得

惯性力忽略不计，式(18)简化为

对式(19)变形、积分后，可得板状管中自发渗吸的隐式解析解tf为

令M=(ρo-ρw)w2gsinα，N=2wσcosθ-ρoLw2gsinα。当x=L时，油水界面到达样品顶端，则单一板状管中渗吸结束的时间tef为

当dt足够小时，误差可以忽略，对式(19)进行差分并积分后，可得板状管中自发渗吸的数值解为

第i根单一板状管的渗吸量(产油量)qfi为

2 多孔介质渗吸数学模型

基于图1岩心尺度渗吸物理模型，连通孔隙空间看作一束曲折的毛细管束，连通微裂缝看作曲折的板状管束，这些曲折毛细管束和板状管束的数目与大小均遵循分形标度规律。基于分形理论[12-15]，对上述渗吸模型进行尺度升级。

2.1 毛细管束模型

根据分形理论，孔隙数与孔径之间的关系为

式中：Np为毛细管半径大于r的孔隙总数；ξ为半径大于r的孔隙半径，m；rmax为最大孔隙半径；Dfp为毛细管的分形维数。

将单元体中岩石样品中半径大于r的毛细管数量按比例放大，得出岩石样品的毛细管数。由YU等[17]推导的表面孔隙度、孔隙面积和横截面的关系，可得岩心样品毛细管数：

式中：d为孔隙直径，m；φp为孔隙所占的孔隙度；τp为毛细管迂曲度[18]，由计算得到，因此毛管长度L=τpH。对式(25)进行微分，计算半径从r到r+dr范围的毛细管数增量：

毛管束中的总流量表达式为[19]

2.2 板状管束模型

板状管开度分布可用分布函数f(w)描述，总管数Nf与分布函数的关系为[20]

根据分形理论，板状管数与开度的关系为

式中：Nf( ≥w)为开度大于w的板状管数量；wmax为最大开度，m；Dff为板状管分形维数。板状管数增量n(w)为

假设所有板状管的迂曲度τp都相同，符合如下关系：

式中：LB为板状管长度，m；H为岩心长度，m。α为板状管的倾斜角。则板状管空间的总体积Vf为

式中：γ为板状管宽度与开度之比，wmin为最小开度，m。同时，微裂缝空间体积占岩心流动空间的关系如下：

式中：D为岩心直径，m。

则由上述公式计算板状管总数Nf为

板状管的概率密度函数为

则板状管总数为

根据分形理论和上述公式，对单条板状管产量进行积分，可以得到岩心尺度微裂缝中自发渗吸的总产油量随孔径分形分布的表达式：

2.3 裂缝性多孔介质渗吸模型

岩心尺度多孔介质包括孔隙和微裂缝2 部分，毛细管和板状管中的自发渗吸存在差别，如图1所示。由上述章节中的公式，可以计算岩心尺度自发渗吸的总产油量，其表达式为

式(38)是考虑孔隙和微裂缝自发渗吸的最终半解析解。图4所示为计算第i个时间步的总渗吸量Q的计算程序，图中下标p 和f 分别表示孔隙和裂缝。利用MATLAB 实现了这个过程。计算孔径范围从rmin到rmax的毛细管数，以及不同开度板状管的概率密度函数。对每个时间步的界面位移距离x进行验证，与岩心长度H进行比较，计算毛管、板状管中油水界面达到岩样顶端时的渗吸量。然后，利用式(27)和式(37)计算特定时间步毛细管束、板状管束的渗吸量。最后，通过将所有渗吸量相加，计算第i个时间步的总渗吸量。根据上述流程计算岩心尺度渗吸量与时间的关系，并对影响因素进行分析。

图4 裂缝性多孔介质自发渗吸总产油量计算流程图Fig.4 Flow chart for calculating total oil production by SI of fractured porous media

3 实例应用与敏感性分析

本节将对上述数学模型进行验证与应用，参照研究区的实际储层参数，开展不同因素对渗吸驱油速度和采出程度的影响研究。利用标度律对岩心尺度渗吸进行尺度升级后，可以结合油藏数值模拟方法提出合理的压后闷井时间建议，对页岩油藏高效开发提供参考。

3.1 岩心尺度渗吸模型应用

页岩油藏储集空间分为孔隙和微裂缝，孔隙以微孔、小孔为主，孔径大都小于1 μm，构成纳米级孔隙系统[21]。根据研究区块的岩心样品，微裂缝较为发育，主要为页理缝和构造微裂缝，成为页岩主要的储集空间和渗流通道[22]。利用扫描电镜对大港油田页岩油藏研究区块的样品进行扫描，页岩储层基质孔隙和微裂缝如图5所示，联合核磁共振技术得到页岩孔径分布特征，其孔隙尺寸和微裂缝开度范围如表1所示。

图5 页岩扫描电镜图Fig.5 SEM of shale

用于模拟和室内实验的岩心与流体参数如表1所示，计算闷井时间与渗吸产油量的关系，开展不同因素对闷井渗吸的影响研究。

表1 岩心与流体物性参数Table 1 Physical parameters of core and fluid

毛细管和板状管的分布对岩心尺度多孔介质中的渗吸速度影响很大。本文利用分形分布理论，进行流体流动空间的敏感性分析。在对岩石流动空间大小进行测定时，常采用压汞法和毛管力曲线研究其微观结构。利用分形定律和计盒法，对像素N(r)所占盒数进行测算：

根据现有的模型[23]，N(r)表达式如下：

式中：VHg为汞侵入岩心的体积。由式(40)可得

根据式(2)，结合式(41)可得

汞饱和度SHg计算公式如下：

由式(42)和式(43)可得

式中：e为常数，结合汞饱和度和毛管对数图可以得到分形维数Df。分形维数可以准确描述储层孔隙结构的不规则性，反映了毛细管和板状管的数量分布。分形维数越大，孔隙半径和微裂缝开度越小，岩心的流通通道越复杂。

根据表1数据和计算流程图(图4)，对不同分形维数模型进行了模拟，图6所示为具有不同分形分布的岩心尺度多孔介质中自发渗吸随时间的产油量。由图6可知：在103s 之前不同分形维数多孔介质的渗吸速度差异不明显，并且速度都很小。在103s 之后，毛细管分形维数和板状管分形维数从1.2增加到1.8，半径小的毛细管和开度小的板状管比例增加，图中3条曲线的斜率依次变小，则页岩储层岩心尺度的渗吸速度显著降低。虽然渗吸速度随分形维数增加而降低，但由于岩心尺度的渗吸结束时间由小孔径流动空间决定，在107s 岩心尺度产油量几乎同时达到平衡。

图6 不同分形分布时岩心渗吸产油量随时间的关系Fig.6 Relationship between oil production of core and time by imbibition with different fractal distributions

图7所示为不同油水黏度比时岩心渗吸产油量随时间的关系。由图7可知：随油水黏度比的增加，渗吸速度降低。图8所示为不同界面张力时岩心渗吸产油量随时间的关系，由图8可知：多孔介质中的渗吸速度随界面张力的增大而增大。本文毛管束和板状管束模型采用分形分布，通过模拟可知油水黏度比和界面张力对岩心尺度渗吸结束时间的影响不明显。

图7 不同油水黏度比时岩心渗吸产油量随时间的关系Fig.7 Relationship between oil production of core and time by SI with different oil-water viscosity ratios

图8 不同界面张力(TIF)时岩心渗吸产油量随时间的关系Fig.8 Relationship between oil production of core and time by SI with different interfacial tension(TIF)

3.2 渗吸实验与模型验证

为了验证裂缝性页岩油藏多孔介质中累计产油量计算式(38)的有效性，利用真实岩心对理论预测进行了检验。针对典型岩心样品开展了渗吸实验，如图9所示，利用核磁共振仪定量表征了孔隙和微裂缝中的渗吸情况，识别了不同渗吸阶段渗流特征。实验与模拟数据如表1所示，其中水相质量分数为8%的KCl水溶液，油相为正十二烷，在室温、大气压下进行的自发渗吸室内实验。同时，根据图4流程图和表1实验参数，进行孔隙和微裂缝分形分布、流动空间内油水界面运移及产油量的计算和表征。

图9 自发渗吸实验示意图与核磁共振仪Fig.9 Schematic diagram of SI experiment and NMR

基于核磁共振T2谱结果，得到了自发渗吸实验岩心产油量随时间的变化关系，如图10所示。同时，图中比较了本文多孔介质渗吸半解析解表达式计算值与渗吸实验数据，由图10可见：渗吸实验观察和理论计算的渗吸产油量与时间几乎匹配一致。

图10 多孔介质渗吸半解析解与渗吸实验结果对比图Fig.10 Comparison of semi-analytical solution of SI and SI experimental results in porous media

2种方法均显示：在渗吸前期产油速度相对较高，毛细管和板状管中的自吸作用是渗吸驱油的最主要作用，尤其是半径小的毛管和开度小的微裂缝对渗吸驱油的速度起重要作用。渗吸中、晚期的产油速度逐渐变小，水相前缘推进到样品端面后驱油速度趋于0。

因此，对于微裂缝发育的页岩油藏，其流动空间相对较小，孔喉连通性一般，闷井过程中的渗吸驱油具有一定的效果，同时较大的微米级孔隙和微裂缝也对渗吸驱油起到了很大的作用。渗吸驱油主要作用力为毛管力，其驱油效率有限，需要通过增加渗吸时间、调整储层压力场等手段提高原油采收率。

4 结论

1) 建立了同时考虑孔隙和微裂缝分形分布的多孔介质渗吸模型，更精确地表征了闷井渗吸产油量随时间的关系。

2) 利用半解析解定量评价了孔隙分布、微裂缝分布和油水黏度等因素对渗吸规律的影响。随油水黏度比减小，界面张力增大，多孔介质中的渗吸速度增大。由数值分析和实验可知自发渗吸的采出程度为17%左右。在页岩油藏压后闷井过程中，增大储层的亲水性，保持一定的油水黏度比和界面张力可以有效发挥渗吸驱油潜力。

3) 渗吸早期产油速度较快，中、后期产油速度逐渐降低，毛管力驱动的自吸采油具有一定的局限性，可通过增加渗吸时间、调整储层压力场、添加表面活性剂等手段控制渗吸效率，达到提高采收率的目的。