一种基于数字成像技术计算致密岩心相对渗透率的方法与流程

本发明涉及致密岩石微观物理性质分析领域，特别是一种基于数字成像技术计算致密岩心相对渗透率的方法。

背景技术：

致密性储层是目前非常规油气资源开发的研究重点，具有广阔的开发前景。油气在低渗致密储层中的流动常伴有非达西渗流现象，为了能有效表征这种低速非达西流动效应，引入了启动压力梯度。由于较低的渗透性和启动压力梯度的存在，使得相对渗透率曲线的测量耗时长准确性差。因此，如何建立起一套考虑多孔介质低渗透特性的相对渗透率计算方法便显得十分关键。

针对致密性储层而言，获取渗透率曲线的方法主要有两种途径：一是开展驱替实验，二是理论计算。由于致密性储层低渗透的特性，使得利用驱替实验来获得相对渗透率曲线的方法耗时长，达到流量稳定所需要的时间较长，且对计量仪器的精度要求较高。对于理论计算方法，又可分为稳态驱替法和动态驱替法两种，他们分别基于稳态达西定律和非稳态jbn理论进行计算。无论是稳态法还是动态法均没有考虑岩石内部结构、孔隙分布等岩石形态特征对相对渗透率曲线的影响，其计算得到的相渗曲线往往与实测曲线相差较大。因此有必要考虑启动压力梯度和毛管力的影响，并结合岩石的孔隙结构特征，找到一种适合于低渗致密岩石的相渗计算方法。

技术实现要素：

为解决现有技术中存在的问题，本发明提供了一种基于数字成像技术计算致密岩心相对渗透率的方法，解决了的问题。

本发明采用的技术方案是，一种基于数字成像技术计算致密岩心相对渗透率的方法，包括以下步骤：

s1：制备满足分辨率要求的致密岩心小柱样；

s2：利用microct-400进行扫描并建立数字岩心；

s3：根据数字岩心，统计分析反映岩石孔隙结构和形状特征的参数；

s4：使用3d图像分形盒维数算法，计算迂曲度分形维数dt和孔隙度分形维数df；

s5：利用标签分析统计最大孔隙等效直径λmax及最小孔隙等效直径λmin；

s6：模拟单根分形毛管内的水驱油过程，求解临界毛管直径λcr，得到驱替时刻t所对应的临界毛管直径；

s7：利用临界毛管直径判断流出端流体类型的方法，计算致密岩心中t时刻所对应的水相体积vw，孔隙体积vp，进一步计算出岩心的含水饱和度sw；

s8：模拟低渗致密岩心中水驱油过程，计算单相流体流量含水饱和度sw不变所对应的单相流体流量qs；

s9：模拟低渗致密岩心中水驱油过程，计算时刻t所对应的流出端油相流量qo以及水相流量qw；

s10：计算驱替时间t所应对的水相相对渗透率krw和油相相对渗透率kro；

s11：改变时间t，判断含水饱和度sw是否不变；若是，则进入s12；若否则返回s6；

s12：绘制出相渗曲线。

优选地，s1的致密岩心直径为5～10mm，长度大于10mm。

优选地，s2包括以下子步骤：

s21：截取目标计算区域；

s22：对目标区域进行滤波和阈值分割，得到二值化图像；

s23：将二值化图像结果重组为三维数字岩心。

优选地，s4的3d图像分形盒维数算法包括以下步骤：

s41：设a为rⁿ空间的任意非空有界子集，对于任意的r＞0，覆盖a所需边长为r的n维立方体的最小数目是nr(a)；

s42：如果存在d，使当r→0时：

nr(a)∝1/r^d

那么称d为a的盒维数，此时存在唯一一个正数k，使得

对上式两边取对数，可得

s43：统计出不同r值时覆盖a分别所需的盒子个数nr(a)，在以lgr为横坐标、以lgnr(a)为纵坐标的对数坐标系中画出(lgr，lgnr(a))，最后通过这些点的拟合线斜率求绝对值，可以得到集合a的分形维数。

优选地，s6包括以下子步骤：

s61：在给定的某一时刻t下，对于毛管等效直径λ≥λcr的单根毛管，只存在单相水流，流出端流体为水；

s62：对于毛管等效直径λ＜λcr的单根等效毛管中存在油水两相流动，流出端流体为油；

s63：将多孔介质中的毛管分成两部分进行分析，λ≥λcr的单根毛管和λ＜λcr的单根毛管。

优选地，s6的临界毛管直径的计算公式为：

式中，jo为油相启动压力梯度；δ为表面张力；θ为接触角；μw为水相粘度；μo为油相粘度；λcr为临界毛细管直径；δp为岩心两端压差；t为时间，s；l为岩心长度；dt为迂曲度分形维数。

优选地，s7中的水相体积vw计算公式为：

孔隙体积vp的计算公式为：

岩心的含水饱和度sw的计算公式为：

式中，dt为迂曲度分形维数；df为孔隙度分形维数；λmax为最大孔隙等效直径；λmin为最小孔隙等效直径，λcr为临界毛细管直径，lt为分形毛管束的实际长度，xt是分形毛管束的直线长度,n表示毛管束的数量。

优选地，s8中单相流体流量qs

式中，dt为迂曲度分形维数；df为孔隙度分形维数；λmax为最大孔隙等效直径；λmin为最小孔隙等效直径，μw为水相粘度，δp为岩心两端压差。

优选地，s9中流出端油相流量qo的计算公式为：

水相流量qw计算式为：

式中，λcr为临界毛细管直径；df为孔隙度分形维数；λcr为临界毛细管直径；μw为水相粘度，jw为水相启动压力梯度。

优选地，s10驱替时间t所应对的水相相对渗透率krw、油相相对渗透率kro，根据单相流达西定律及两相流达西定律可得，

式中，qw为水相流量，qo为流出端油相流量，qs为单相流体流量，μw为水相粘度，μo为油相粘度。

本发明基于数字成像技术计算致密岩心相对渗透率的方法的有益效果如下：

1.该方法考虑了启动压力梯度和毛管压力的影响，体现了致密岩心低渗的特点。

2.该方法基于多孔介质分形毛管数模型进行相对渗透率计算，避免了实验的耗时和仪器精度问题。

3.该发明基于室内ct扫描实验结果，建立对应的真实岩样数字岩心，考虑了岩石内部结构、孔隙分布等岩石形态特征对相对渗透率曲线的影响，并且可进行反复多次的计算模拟，避免了重复物理实验的资源浪费。

附图说明

图1为本发明基于数字成像技术计算致密岩心相对渗透率的方法的提供的致密岩心相对渗透率计算流程图。

图2为本发明基于数字成像技术计算致密岩心相对渗透率的方法实施例得到的数字岩心图。

图3为本发明基于数字成像技术计算致密岩心相对渗透率的方法的计算涉及的单根分形毛管流动示意图。

图4为本发明基于数字成像技术计算致密岩心相对渗透率的方法的计算涉及的岩心水驱油分形毛管流动示意图。

图5为本发明基于数字成像技术计算致密岩心相对渗透率的方法的实施例1得到的相渗曲线与实测数据的对比图。

图6为本发明基于数字成像技术计算致密岩心相对渗透率的方法的实施例2得到的不同分形维数对应的相渗曲线对比图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种基于数字成像技术计算致密岩心相对渗透率的方法，包括以下步骤：

s1：制备满足分辨率要求的致密岩心小柱样；

s2：利用microct-400进行扫描并建立数字岩心，如图2所示；

s3：根据数字岩心，统计分析反映岩石孔隙结构和形状特征的参数；

s4：使用3d图像分形盒维数算法，计算迂曲度分形维数dt和孔隙度分形维数df；

s5：利用标签分析统计最大孔隙等效直径λmax及最小孔隙等效直径λmin，如图4所示；

s6：模拟单根分形毛管内的水驱油过程，求解临界毛管直径λcr，得到驱替时刻t所对应的临界毛管直径；

s7：利用临界毛管直径判断流出端流体类型的方法，计算致密岩心中t时刻所对应的水相体积vw，孔隙体积vp，进一步计算出岩心的含水饱和度sw；

s8：模拟低渗致密岩心中水驱油过程，计算单相流体流量含水饱和度sw不变所对应的单相流体流量qs；

s9：模拟低渗致密岩心中水驱油过程，计算时刻t所对应的流出端油相流量qo以及水相流量qw；

s10：计算驱替时间t所应对的水相相对渗透率krw和油相相对渗透率kro；

s11：改变时间t，判断含水饱和度sw是否不变；若是，则进入s12；若否则返回s6；

s12：绘制出相渗曲线。

本实施方案的s1的致密岩心直径为5～10mm，长度大于10mm。

本实施方案的s2包括以下子步骤：

s21：截取目标计算区域；

s22：对目标区域进行滤波和阈值分割，得到二值化图像；

s23：将二值化图像结果重组为三维数字岩心。

本实施方案的s4的3d图像分形盒维数算法包括以下步骤：

s41：设a为rⁿ空间的任意非空有界子集，对于任意的r＞0，覆盖a所需边长为r的n维立方体的最小数目是nr(a)；

s42：如果存在d，使当r→0时：

nr(a)∝1/r^d

那么称d为a的盒维数，此时存在唯一一个正数k，使得

对上式两边取对数，可得

s43：统计出不同r值时覆盖a分别所需的盒子个数nr(a)，在以lgr为横坐标、以lgnr(a)为纵坐标的对数坐标系中画出(lgr，lgnr(a))，最后通过这些点的拟合线斜率求绝对值，可以得到集合a的分形维数。

如图3所示，本实施方案的s6包括以下子步骤：

s61：在给定的某一时刻t下，对于毛管等效直径λ≥λcr的单根毛管，只存在单相水流，流出端流体为水；

s62：对于毛管等效直径λ＜λcr的单根等效毛管中存在油水两相流动，流出端流体为油；

s63：将多孔介质中的毛管分成两部分进行分析，λ≥λcr的单根毛管和λ＜λcr的单根毛管。

本实施方案的s6的临界毛管直径的计算公式为：

式中，jo为油相启动压力梯度；δ为表面张力；θ为接触角；μw为水相粘度；μo为油相粘度；λcr为临界毛细管直径；δp为岩心两端压差；t为时间，s；l为岩心长度；dt为迂曲度分形维数。

本实施方案的s7中的水相体积vw计算公式为：

孔隙体积vp的计算公式为：

岩心的含水饱和度sw的计算公式为：

式中，dt为迂曲度分形维数；df为孔隙度分形维数；λmax为最大孔隙等效直径；λmin为最小孔隙等效直径，λcr为临界毛细管直径，lt为分形毛管束的实际长度，xt是分形毛管束的直线长度,n表示毛管束的数量。

本实施方案的s8中单相流体流量qs

式中，dt为迂曲度分形维数；df为孔隙度分形维数；λmax为最大孔隙等效直径；λmin为最小孔隙等效直径，μw为水相粘度，δp为岩心两端压差。

本实施方案的s9中流出端油相流量qo的计算公式为

水相流量qw计算式为：

式中，λcr为临界毛细管直径；df为孔隙度分形维数；λcr为临界毛细管直径；μw为水相粘度，jw为水相启动压力梯度。

本实施方案的s10驱替时间t所应对的水相相对渗透率krw、油相相对渗透率kro，根据单相流达西定律及两相流达西定律可得，

式中，qw为水相流量，qo为流出端油相流量，qs为单相流体流量，μw为水相粘度，μo为油相粘度。

本实施方案在实施时，实施例一

(1)选取低渗致密砂岩岩样，制备满足分辨率要求的小柱样，尺寸大小为8mm(直径)×20mm，然后利用microct-400进行扫描并建立数字岩心，如图2。

(2)使用3d图像分形盒维数算法计算迂曲度分形维数dt＝1.1及孔隙度分形维数df＝1.3，利用标签分析统计最大孔隙等效直径λmax＝2.6×10^-5m及最小孔隙等效直径λmin＝4.3×10^-8m。

(3)以t＝0作为起始时刻，计算对应的临界毛管直径，判断该时刻岩心内所有毛管流体分布情况。

(4)模拟计算单相流动的流量qs，然后模拟计算t＝0时水驱油流动时的水相流量qw及油相流量qo。

(5)模拟计算t＝0时岩心内的水相体积、孔隙体积及含水饱和度；

(6)计算该含水饱和度对应的相对渗透率kro及krw。

(7)重复事实例一的步骤(3)～(6)，计算下一个时间步t＝1的各项参数，然后t+1，直至计算的岩心含水饱和度不变，绘制出相渗曲线，如图5。

将计算所得的相渗曲线与实测数据进行对比，可见计算可靠度高，上述实施例有力地证明了本发明的有益效果。

实施例二

(1)选取不同的孔隙度分形维数df＝1.3、1.4、1.5，相同的迂曲度分形维数dt、最大孔隙等效直径λmax及最小孔隙等效直径λmin，作为岩心孔隙结构参数的特征值，用以计算不同孔隙度分形维数对相对渗透率曲线的影响。

(2)以t＝0作为起始时刻，计算对应的临界毛管直径，判断该时刻岩心内所有毛管流体分布情况。

(4)模拟计算单相流动的流量qs，然后模拟计算t＝0时水驱油流动时的水相流量qw及油相流量qo。

(5)模拟计算t＝0时岩心内的水相体积、孔隙体积及含水饱和度；

(6)计算该含水饱和度对应的相对渗透率kro及krw。

(7)重复实施例二的步骤(3)～(6)，计算下一个时间步t＝1的各项参数，然后t+1，直至计算的岩心含水饱和度不变，绘制出相渗曲线，如图6。