基于格子玻尔兹曼方法的岩心切片渗流通道研究

2020-12-29刘博伟牟松茹李彦来张文童

石油化工高等学校学报 2020年6期

刘博伟，牟松茹，李彦来，张墨，张文童

（中海石油（中国）有限公司天津分公司渤海石油研究院，天津300459）

数字岩心技术是一种依托高精度成像设备在岩心孔隙级层面建模，继而分析储层物性参数的技术，目前已成为一种较为成熟、可用于油气田开发工程的分析手段。相对于传统室内岩心实验，数字岩心以其实验周期短、实验可重复、微观因素可控、节约岩心等优势被广泛关注及应用。数字岩心扫描分辨率已达到纳米级，将油气流动研究推向微观孔隙流动研究的大门[1-5]。

数字岩心孔隙流动具有异常复杂的孔隙空间，孔隙和骨架的边界极不规则，直接求解Navier-Stokes方程不太符合实际，有限元、有限体积、有限差分等方法也难以处理异常复杂边界[6-7]。格子玻尔兹曼方法（Lattice Boltzmann Method，LBM），作为介观尺度计算流体力学方法，具有介于微观分子动力学模型和宏观连续模型的介观模型特点，因其流体相互作用易描述、复杂边界易设置、并行计算容易等优势，逐渐被应用到多孔介质复杂空间流动模拟[8-13]。

本文以渤海S油田东营组储层数字岩心切片为研究对象，通过图像处理完成孔隙网络建模，利用格子玻尔兹曼方法开展岩心切片二维孔隙流动模拟及速度场表征，在此基础上设计实验，分析生产压差、孔喉优化对渗流通道的影响，探寻合适方法提高微观孔隙渗流波及程度。

1 数字岩心切片流动模拟

实验采用渤海S油田东营组某砂岩岩心，岩心储层深度为1 379.0～1 379.1 m，岩心CT扫描结果如图1所示，其中黑色部分为孔隙，灰色部分为骨架。图1（a）为全直径岩心扫描结果，对应扫描分辨率为 100 μm；图 1（b）为岩心柱塞扫描结果，对应扫描分辨率为 13 μm；图 1（c）为柱塞子样扫描结果，对应扫描分辨率为2 μm。

图1 渤海S油田某储层数字岩心扫描结果Fig.1 Digital core scanning results of a reservoir in Bohai S oilfield

根据岩心扫描切片可知，全直径岩心由于扫描精度很低，无法有效识别岩心孔隙及裂缝；岩心柱塞扫描分辨率中等，岩心骨架边缘灰度相近，识别存在误差，影响岩心孔隙分析；柱塞子样扫描分辨率较高，能够有效识别骨架边缘、岩心孔隙。因此，本文基于柱塞子样扫描结果，随机选取岩心某一切片开展孔隙网络建模及切片流动模拟（见图2），分析切片渗流情况。

图2 数字岩心某二维切片Fig.2 Digital two-dimensional slice

1.1 岩心切片图像处理

在数字岩心分层CT扫描切片过程中，受扫描设备影响不可避免地会产生一定程度的噪声，干扰骨架边界识别，因此需要开展灰度图像降噪处理。常用降噪处理方式有均值滤波和中值滤波。由于中值滤波属非线性降噪技术，能够相对有效地保护边缘信息，在数字岩心灰度图像降噪处理的同时兼顾保护骨架边界的不规则形状。因此，选用中值滤波技术对数字岩心扫描图像进行噪声过滤，降噪效果见图3。

图3 切片局部原始及过滤后灰度图像Fig.3 Local original and filtered grayscale image of a core plunger

数字岩心CT扫描结果为灰度图像堆栈，而用于流动模拟的孔隙网络模型仅关注骨架、孔隙两种状态，因此需要开展阈值分割，将灰度图像二值化，用1代表骨架，0代表孔隙。采用单阈值处理方式，数字岩心数据体设定单一阈值，灰度值大于阈值视为骨架，小于阈值视为孔隙。选用试凑法以储层孔隙度为依据进行阈值选取，当阈值取76时，数字岩心孔隙度与储层实际值一致，以此为阈值进行图像二值化，二值化处理效果见图4。

图4 数字岩心切片图像二值化处理Fig.4 Digital core slice image binarization

1.2 格子玻尔兹曼方法模型

基于Y.H.Qian等[14]所提出的DdQm系列格子玻尔兹曼方法模型开展二维切片流动模拟，本文仅选用D2Q9系列模型，通过C++语言编程实现模型计算。

格子玻尔兹曼模型主要由三部分组成，即格子（离散速度模型）、平衡态分布函数、分布函数的演化方程[15-16]。

数字岩心切片为二维二值化图像，因此离散速度模型选用D2Q9模型（即适用范围为二维空间，模型包含9个离散速度方向，见图5），离散速度配置见式（1）。

图5 D2Q9离散速度模型Fig.5 D2Q9 discrete velocity model

其中，α=0,1,2,…,8代表D2Q9模型离散速度方向；eα代表对应α方向上的速度；c代表单位格子速度，即格子步长与时间步长之比。

平衡态分布函数选用Y.H.Qian等[14]提出的LBGK模型（由于本文所研究的切片渗流仅局限于单相流体流动，因此不考虑附加作用力项），平衡态分布函数见式（2）、（3）。

平衡态分布函数对应的演化方程见式（4）。

其中，r代表流体粒子的空间位置；t代表时间；δt代表时间步长。

根据LBGK模型理论，宏观密度、速度定义见式（5）。

其中，ρ′和u′分别表示模型的宏观密度、速度。

数字岩心切片流动模拟对应的孔隙网络模型呈矩形（见图6），入口端/出口端选用压力边界，顶面/底面及孔隙骨架边界选用标准反弹边界。标准反弹边界，指在某微粒格点上，微粒按照运动轨迹运动至边界格点，其后微粒发生反弹，反弹后沿原运动轨迹相反的方向运动；标准反弹边界常用于描述无滑移壁面反弹，即沿壁面方向无速度梯度的反弹。压力边界条件，参考Q.S.Zou等[17]根据非平衡反弹思想提出的动力学压力边界格式假设；在孔隙网络模型流入端压力值已知，即 Pin=cs2ρin（Pin为流入面压力，ρin为流入微粒密度）已知，而流入端速度未知,边界处沿壁面方向速度分量为零，水平速度分量ux未知，结合离散速度模型及平衡态分布函数与宏观变量的关系，联立方程组即可求解边界速度。

图6 孔隙网络模型示意Fig.6 Schematic diagram of pore network model

1.3 切片渗流过程表征

基于格子玻尔兹曼方法切片流动模拟结果，计算孔隙内流体速度场分布情况，通过对比渗流速度大小区分判断孔隙优势渗流通道。图7为切片孔隙内流体速度场随时间变化情况。当速度场不再随时间发生变化时，达到稳态。本文优势通道研究均是在渗流达到稳态条件下开展的，选取相邻时间步长速度变化量小于10-3为稳态判断条件。

2 岩心切片优化渗流通道分析

基于格子玻尔兹曼方法，结合达西渗流理论，开展不同生产压差、不同孔喉尺寸条件下流体稳态速度场分布研究，分析二维切片优势渗流通道的影响因素，开展优化孔隙渗流方法探索。

数字岩心二维切片宏观尺寸为1.00 mm×0.56 mm，格子尺寸为500×280；实验模拟流体为气体，宏观密度为1.293 kg/m3。

2.1 生产压差对渗流影响

设置入口端/出口端生产压差分别为0.005、0.010、0.015、0.020、0.025 MPa，其他条件不变，可获得对应稳态速度场分布，图8为不同压差条件下渗流通道表征图。

根据图8可知，伴随生产压差逐渐增大，切片孔隙内的渗流速度逐渐增大，流体的渗流通道也趋于稳定。图8（a）中，由于生产压差较小，渗流速度较小，但能够辨别出渗流通道为主通道；图8（b-e）中，渗流通道仍以主通道为主，但随着压差进一步增大，次级通道也有较弱的渗流。

图7 岩心二维切片速度场Fig.7 Core two-dimensional slice velocity field diagram

图8 不同压差条件下渗流通道图Fig.8 Seepage channel diagram under different pressure difference conditions

根据模拟结果可知，仅通过放大生产压差，次级通道渗流速度较小，无法高效驱替局部微观剩余油。

2.2 孔喉优化对微观渗流影响

利用数字岩心微观因素可控的优势，结合二维切片优势渗流通道展布形态，类比现有的油田增产技术（储层酸化压裂，即扩大岩心孔缝；注采井网调驱，即封堵孔喉改变流场），有针对性地改造二维切片局部骨架形态，测试其对孔隙微观渗流的影响。

假设能够实现储层岩石微观骨架孔喉封堵，以此为基础开展二维切片骨架改造。将原二维切片主通道入口进行孔喉封堵，结果如图9所示。根据稳态渗流场分布可知，孔喉封堵能够有效改变切片渗流通道分布，难渗流区域开始形成有效渗流，达到增大微观孔隙波及程度的目的。然而，次级通道入口孔径较小，封堵优势通道也造成渗流速度变小，不利于储层渗流。

图9 封堵局部孔喉条件下渗流通道Fig.9 Seepage channel diagram under local pore throat conditions

针对微观孔喉封堵后造成的渗流速度变小的问题，提出微观压裂酸化的假设，以此进行切片骨架孔喉扩径改造，增大孔隙渗流速度。在封堵主通道入口的基础上，将局部区域孔喉扩径，增大次级通道入口的半径，结果如图10所示。根据稳态渗流场分布可知，难渗流区域开始形成高效渗流，切片微观波及程度进一步扩大，达到扩大微观储层改造的目的。

基于二维切片骨架改造的假设，仅针对图10中次级通道入口进行改造，开展不同孔径大小对孔隙渗流速度的影响。模拟结果如图11所示，当次级通道入口孔径相对较小时，其对渗流速度的制约作用越大；当孔径扩大到一定程度时，渗流速度增幅逐渐减小。