基于数字岩心的低渗透率储层微观渗流和电传导数值模拟

王平全, 陶鹏, 刘建仪, 黄丽莎

(1.西南石油大学油气藏地质及开发工程国家重点实验室, 四川 成都 610500; 2.长江大学资源与环境学院, 湖北 武汉 430100)

0 引 言

合理进行页岩气勘探开发,必须研究其孔隙微观渗流特征。孔隙微观渗流机理研究是在微米尺度上研究流体的各种流动规律以及与宏观性质之间的联系。随着计算机技术的发展,目前利用数字岩心三维重构研究渗流机理已经成为一个热门的研究方向[1-3],实现孔隙级别上的数据可视化,将有助于对各种流动机理的认识。由于CT技术具有无损检测特性,用它获取岩心二维切片及三维立体图像是当前建立数字岩心三维重构可视化最精确最直接的办法。本文利用Avizo数据可视化软件可以形象生动地展现所研究问题的实质,同时与COMSOL有限元软件对接,能有效解决当前低渗透率储层微观渗流问题。在孔隙纳米级渗流机理研究中利用该软件成功地实现了各种数据体的可视化及渗流研究。

1 数字岩心三维重构的可视化

在微米级渗流机理研究中,精细表征孔隙空间结构是亟待解决的重大问题。目前主要有2大类方法:物理实验法和数值重建法。物理实验法就是借助于高精度的光学仪器获取岩心的平面图像,序列成像法、聚焦扫描法和CT扫描法是应用较为广泛的物理实验方法。数值重建法就是借助岩心的二维切片,利用各种算法获取岩心的孔隙结构特征进行数字岩心三维重构,按所采用算法的不同,数值重建方法又可以分为2类:随机方法和过程模拟法,其中随机法主要包括高斯随机场法、模拟退火法、顺序指示模拟法、多点地质统计学方法和马尔科夫链方法等[4-7]。数值重建法具有实验成本低建模速度快的优点,但这类方法重建数字岩心的连通性较差,反映的岩心孔隙结构与真实多孔介质有一定差距。

利用Avizo对CT扫描法所重建的数字岩心进行了可视化。从简单的可视化和分析测量到先进复杂的图像处理、量化和骨架化,Avizo都为2D和3D可视化、基于有限元分析的3D建模和物理属性的计算等提供了一个全面、多模态的数字实验室。

1.1 Avizo中的图像处理

一个基本的图像处理流程包括载入数据、滤波、二值化处理、图像分割以及渗流方向上的孔隙连通性判断,其中最重要的是滤波和二值化分割。

(2) 图像的二值化处理。二值化处理是一个重要的前期图像处理过程,是后续图像处理技术的基础,本文所用到的岩心切片图像包括孔隙喉道、岩石骨架以及实验过程中的噪声,从多值的数字图像中直接提取出孔隙和骨架信息,最常用的方法是设定一个阈值T,用T值把图像的数据分成2部分:大于T的像素群和小于T的像素群,这种方法称为图像的二值化。具体方法是将图像上的像素点的灰度值设置为0或255,也就是将整个图像呈现出明显的只有黑和白的视觉效果。基于二值化的阈值分割是一种最常用的并行区域技术,它是图像分割中应用数量最多的一类。图像的二值化处理已经确定了阈值T,将阈值与各像素点的灰度值进行比对,二值化分割的结果直接将图像进行分割,一般在Avizo中添加Interactive Thresholding命令模块进行图像分割。

1.2 数值模拟结果

经过一系列图像处理流程,得到了孔隙度、孔喉大小和方位、孔隙空间面积、骨架体积等的定量研究结果,达到对岩心微观孔隙结构的定量描述[8-9]。由于得到的数据量巨大,通过统计每个孔隙面积的处于某个范围内的个数描述孔隙结构特征,利用Avizo软件得到了岩心三维重构和孔隙网络模型结果示意图(见图1)。

图1 岩心微观孔隙结构定量表征

2 Avizo与COMSOL数据交互

数字岩心孔隙空间几何模型三维重构是为了进行三维数据的可视化,也是为微观渗流模拟等后期应用做准备,因此当数字岩心建立后选择合适的数据接口对后续的研究有着非常重要的意义。一般地,作为交互对接的数据接口有2类,一类是OFF文件类型;另一类是STL文件类型。Object File Format(OFF)文件是一种能够通过描述物体表面的多边形表示几何模型的文件类型,OFF文件结构简单,需要的存储量小但通用性差。STereo Lithography(STL)文件具有结构简单、通用性强、易于实现等优点,是数字岩心三维重构比较常用的文件格式之一。但在实际应用中,STL文件包含的三角网格数量较多,文件存在信息冗余现象,因此在使用STL文件类型时,必须经过平滑、简化、修复表面和增强网格质量等一系列流程,本文主要采用STL数据接口实现Avizo数字岩心可视化软件和COMSOL有限元模拟软件的交互对接。为得到高质量的网格并运用于后期有限元模拟中,需要对Avizo创建的数字岩心表面网格质量进行改进,具体的做法包括数据体表面简化、三角形表面修复和三角网格优化。

2.1 数据体表面简化

当孔隙结构比较复杂、数字岩心尺寸较大时,直接生成STL文件会占据计算机很大储存空间和计算资源,需要对数据体表面进行简化。在平面区域中生成大三角形,同时在高曲率区域生成较小的三角形,这样保留细节导出轻量级的三角形表面文件。

一般地,太强烈的简化可能导致三角形网格相交从而产生四面体网格,此时应该再次执行网格质量检测,进一步进行表面修复和网格优化,排除可能的表面异常从而增强三角网格质量。

2.2 三角形表面修复和三角网格优化

低渗透砂岩数字岩心的孔隙空间异常复杂,进行三维重构提取表面和划分三角网格时存在困难,模型完成网格划分后,由于网格质量很差,在后续导入有限元软件进行数值模拟时会产生诸多错误,无法进行渗流模拟。采用Avizo可视化软件所带的表面编辑器模块对数字岩心孔隙空间模型的表面细节部分进行人机交互式修复操作,其中包括修复缺口、尖角,消除交叉、缝隙、倒角、共面、小孔、重合边等,最终生成质量较好的三角网格(见图2)。

图2 孔隙几何模型表面细节修复

完成上述人机交互操作后,通过Avizo网格测试的孔隙几何空间模型可转化为COMSOL有限元软件所能识别的STL文件格式,架起了数字岩心几何模型到微观渗流模拟的桥梁,实现数据交互对接。

3 Avizo与COMSOL对接在模拟微观渗流研究中的应用

利用可视化软件Avizo对岩心进行图像分割并创建孔隙空间表面,然后将其用于有限元多物理场软件COMSOL中。在CT扫描实验中使用的岩样为长度0.45 cm、直径0.55 cm的圆柱体,扫描的体素大小为5.69 μm。

图3表示Avizo中孔隙空间三维重构图像区域,不同簇的连通孔隙以不同的颜色显示,样品中只有大的绿色区域相互连通,允许流体和电流流过模型域(假定基质为非导电矿物),对于岩心样品,仅选择该区域用于模拟。

图3 Avizo中孔隙空间3D图像

图4表示COMSOL中孔隙空间三维网格,实验中考虑垂直于立方体岩心的边界(即模型域)方向上的流动。当孔隙表面与其法线平行于给定模拟的流动方向并与边界相交时,认为是入口或出口区域的一部分。

图4 COMSOL中的孔隙网格的3D图像

将STL文件导入到有限元软件COMSOL中,需要STL几何修正校正某些模拟的拓扑。对于这些模拟,需要完成COMSOL中的面分割,包括调整孔隙表面中的最大角度、最大面邻近角度、最小相对面积和小面的移除。在一些特殊的情况下,有必要返回到Avizo中进行更多的手动调节STL文件,才能成功地融合在COMSOL中,完成手动操作后,COMSOL能够生成一个没有内部几何误差的自由四面体3D网格。

3.1 理论基础

渗流模拟是基于COMSOL层流模块中稳态不可压缩Navier-Stokes方程,包括力学平衡方程[10-11]

ρu·u=-p+·η[u+(u)T]

(1)

还涉及到质量平衡方程

·u=0

(2)

式中,ρ和η分别为流体密度和流体动态黏度;p和u分别为流体的压力和速度。爬流通常发生在自然压力梯度下的多孔岩石,对这样的流动而言,流体运动方程中的动量项很小,即式(1)中的左边可以忽略不计,力学平衡主要受压力梯度项和黏性应力影响,即式(1)中的右边两项。

对于电导率模拟,主要使用电流模块,在孔隙空间求解岩石电导率ν的方程是电势的拉普拉斯方程

2ν=0

(3)

假设岩石的电导率和孔隙流体不变,通过岩石基质的电传导忽略不计,电流密度与电势相关

J=-σV

(4)

假定流动模拟入口处的边界条件是恒定压力且没有黏性正应力,即p=1 Pa,μη[u+(u)T]n=0,n表示入口处的一个单位正向量。对于电导率模拟,入口处的边界条件:ν=1 V;出口处,设定p=0,μη[u+(u)T]n=0。出口处的边界条件:ν=0,假设在孔隙的其他面上不存在滑移现象,即u=0且孔隙的其他面是电绝缘的,ν·n=0。

模拟过程中使用的参数值:ρ=1 kg/m3,η=1 Pa·s,σ=1 S/m,通过对入口和出口处的速度和电流密度法向分量的表面积分计算总质量流量Q和电流I,入口和出口流量应该相同,并且差值用做模拟误差的测量。

对于爬流状态下流体流动模拟,质量流量与压力梯度成线性关系且流速与压降成比例,使用达西定律计算样品的渗透率K

(5)

式中,L和A是岩石样品在流动方向上的长度和垂直于压力梯度的岩心横截面积;Δp是岩石样品上的压降。类似地,地层因素F由欧姆定律计算

(6)

Kozeny-Carman方程认为渗透率的变化与孔隙度的变化相关联。在低孔隙度低渗透率岩石样品中,这种关系近似视为幂律接近3的幂。

一般运用Archie公式(Archie,1942)将岩层的形成因子与其孔隙度相关联,具体方程为

F=aφ-m

(7)

常数a通常接近1,而m被称为胶结指数,其对于碳酸盐岩(Doveton,1999)和砂岩(Dvorkin等,2011)通常接近2。

3.2 渗透率和地层因素的模拟结果

图5显示了岩心在y方向施加压力梯度的流体流动的模拟结果,流线平行于流体流动方向,颜色表示速度大小。一些孔隙空间不包含流线,表明这些孔隙对渗透率没有显著影响;还可以看出,流速最大的区域,孔隙空间变得收缩,由于孔壁的阻力,穿过孔壁垂直于流动方向的流速具有明显的梯度。

图5 岩石样品中流体流动模拟结果

图6 岩石样品中电传导模拟结果

在图6中,得到了同一块岩心电导率的模拟结果,其中流线平行于电流方向,颜色指示电流密度的大小。电压施加在y方向上,流动路径可以看做类似于流体流动的情况。然而,在垂直于流体流动方向上的电流密度没有梯度,这是基于流体流动和电流流动边界条件的不同特性所得到的正常电流密度。对于流体流动的情况,孔壁在流体上施加阻力,其与流体黏度共同导致流过孔隙的流速梯度。对于电流流动的情况,孔壁处的切向电流密度不受边界条件约束。

4 结 论

(1) 基于Avizo的图像处理过程形象地揭示出数字岩心的基本问题:微观表征岩心的孔渗饱特性以及流体渗流特征,孔隙级渗流机理研究是在孔隙尺度上进行多相渗流机理研究。借助科学的可视化软件,可以更为形象地反映出研究问题的本质。

(2) 通过Avizo软件与有限元模拟软件COMSOL、Ansys等的结合能有效进行数字岩心的研究,为开展油藏精细描述和油气藏微观渗流机理研究奠定了良好的基础。

参考文献:

[1] BLUNT M J. Flow in Porous Media-pore-network Models and Multiphase Flow [J]. Current Opinion in Colloids and Interface Science, 2001, 6(3): 197-207.

[2] 陶军, 姚军, 赵秀才. 利用IRIS Explorer数据可视化软件进行孔隙级数字岩心可视化研究 [J]. 石油天然气学报, 2006, 28(5): 51-53.

[3] 万金彬, 孙宝佃, 陈守军, 等. 数字岩心技术研究及应用 [J]. 测井技术, 2012, 36(2): 154-159.

[4] 赵秀才. 数字岩心及孔隙网络模型重构方法研究 [D]. 北京: 中国石油大学, 2009.

[5] 姜黎明, 刘宁静, 孙建孟, 等. 利用CT图像与压汞核磁共振构建高精度三维数字岩心 [J]. 测井技术, 2016, 40(4): 404-407.

[6] ØREN P E, BAKKE S. Reconstruction of Berea Sandstone and Pore-Scale Modeling of Wettability Efects [J]. J Pet Sci Eng, 2003, 39: 177-199.

[7] HOU J, ZHANG S K, SUN R Y, et al. Reconstruction of 3D Network Model Through CT Scanning [C]∥Proceedings of the SPE Europe/EAGE Annual Conference and Exhibition. Richardson: Society of Petroleum Engineers, 2007, 1-10.

[8] 陈书荣, 王达健, 张雄飞. 多孔介质孔隙结构的网络模型应用 [J]. 计算机与应用化学, 2001, 18(6): 531-535.

[9] 李菊花, 郑斌. 微观孔隙分形表征新方法及其在页岩储层中的应用 [J]. 天然气工业, 2015, 35(5): 52-59.

[10] 张思勤, 汪志明, 洪凯, 等. 基于格子Boltzmann方法的3D数字岩心渗流特征分析 [J]. 测井技术, 2016, 40(1): 18-22.

[11] 闫国亮, 孙建孟, 刘学锋, 等. 储层岩石微观孔隙结构特征及其对渗透率影响 [J]. 测井技术, 2014, 38(1): 28-32.