非常规储层孔隙网络两相流动模拟研究

王秀坤,刘海成,吴忠维,崔传智

(1.中国石油大学(北京)非常规油气科学技术研究院，北京 102249；2.中国石化胜利油田分公司勘探开发研究院，东营 257015；3.中国石油大学(华东)石油工程学院，青岛 266580)

油气藏多相流动模拟的关键在于表征其毛管力与相对渗透率曲线，传统油气藏主要依靠室内岩心实验测定获取，一般耗时较长，尤其对于非常规储层，如致密、页岩储层，其实验往往难以开展。基于数字岩心孔隙网络流动模拟技术，在孔隙尺度上高效模拟非常规储层内油、气、水的多相流动过程，并计算毛管力与相对渗透率曲线，是行之有效的研究方法之一。

借助微纳米CT(computed tomography)、FIB-SEM(focused ion beam scanning electron microscope)等实验手段，可获得岩石微观孔隙空间的三维数据体，二值化分割后即可获得数字岩心[1]。孔隙尺度流动模拟方法分为直接模拟和孔隙网络模拟，前者直接在数字岩心的孔隙空间的体素点上开展数值模拟计算[2-3]，计算成本高。而后者则一般利用最大球体法[4]、分水岭算法[5]等算法抽提获得表征其孔喉几何形状和拓扑空间的孔隙网络，然后据此开展多相流动模拟计算[6]。通常采用的准静态孔隙网络流动模拟方法只考虑毛细管力的作用，可快速模拟排驱和吮吸过程中油水的三维空间展布变化，继而获取宏观的流动物性参数，如绝对渗透率、相对渗透率曲线和毛细管力曲线等。起初，孔隙网络两相流动模拟法主要是在挪威国家石油公司研究中心的Oren等[7]和英国帝国理工大学Blunt[8]研究基础之上而创建的。Valvatne等[9]于2004年开发了业界公认的一套的两相流动模拟软件，实现了对初次排驱和二次吮吸过程的模拟，可得到相对渗透率、毛细管力曲线和电阻率变化规律等宏观参数，极大推动了孔隙网络流动模拟的研究工作。但 van Dijke等[10]指出，Blunt软件中对油膜坍塌机制的考虑不合理，并利用能量守恒提出了更精确的油膜存在条件模型。Helland等[11]指出了Blunt软件中斜三角形截面传导系数的误差，并提出了更贴近孔喉特征的四角星形截面模型。Ryazanov[12]据此四星形模型，建立了常规储层的孔隙网络两相流动模拟方法，而对非常规储层没有涉及。Petrovskyy等[13]完善了Ryazanov[12]孔隙网络多相流动的模拟计算速度，提出了更高效的非润湿相圈闭检测算法。近年来，中国学者也对孔隙网络流动模拟方法开展了与国际接轨的相关性研究，取得了一定的特色成果[14-16]。

目前，基于孔隙网络模型的非常规储层两相流动模拟研究较少。为此，基于四角星形截面的处理方式，结合高性能图论算法和开源可视化表征软件Paraview[17]，在考虑活塞式驱替、卡断填充、协同式填充、带油膜充填等孔喉微观填充机理基础上，建立了孔隙网络两相流动模拟方法；通过与实验结果对比验证了该方法的有效性和准确性；并运用该孔隙网络两相流动模拟方法对吉木萨尔页岩储层的复杂油水两相流动机理进行了探讨。

1 孔隙网络两相流动的模拟方法

1.1 基本几何关系及排驱过程模拟

当孔隙网络模型中每个孔和喉的截面采用四角星形时(图1)，通过几何关系可得孔或喉的截面积A，其计算公式为

图1 四角星形截面几何关系示意图Fig.1 Schematic of the four-pointed crossing sectional area

A=2r2(1+cotγ)

(1)

式(1)中：A为孔或喉的截面积，m2；r为孔或喉内切圆半径，m；γ为四角星的半角，(°)。

(2)

(3)

式中：G为截面的形状因子，无量纲；P为周长，m。

对于初次排驱过程，即油藏成藏过程，储层首先被水充满，烃源岩生油后通过生烃压力的作用进入储层，油相克服毛细管力驱替水相，最终形成油藏。排驱和后续吮吸过程，均认为岩石为水湿。对于这个过程的孔隙网络模拟，毛细管力逐步增大，喉道起到主要的约束作用，任一喉道入口毛细管力大小可表示为式(4)[12]。当毛细管力超过式(4)时，油相进入此喉道以及与之连接的孔隙。

(4)

每个孔或喉的含水饱和度Sw可通过截面积占比获得，可表示为

(5)

1.2 水湿油藏水驱油过程模拟

对于水湿油藏水驱油过程，即二次吮吸过程，毛细管力逐步降低，水相逐步填充孔和喉。喉道主要有活塞式填充和卡断填充两种微观填充机理，活塞式喉道填充发生时需要其所连接的孔隙有一个已被水相填充，否则只能发生卡断填充。喉道卡断填充的临界毛细管力为[12]

(6)

需要说明的是，当孔隙被水相填充后，与之相连接的喉道也会随之被填充，孔隙制约着喉道活塞式充填的发生，这与排驱过程中喉道制约着孔隙被油相侵入的道理是相对应的。孔隙一般为协同式填充，其临界毛细管力为[9]

(7)

式(7)中：

对于排驱和吮吸的过程，需要考虑润湿滞后现象[18]。排驱过程中采用后退角，而吮吸过程采用前进角，且前进角大于后退角。初次排驱结束之后，标志着油藏的形成和最大的毛细管力。二次吮吸过程发生时，毛细管力逐步降低，此时角隅的水相区域膨胀，接触角从后退角逐步增大，直至变为前进角，含水饱和度缓慢增加，但仅以角隅和水膜的形式存在，对应的含水饱和度计算公式见式(5)。一般在模拟计算时会直接给定前进角和后退角的值，水湿岩心，典型的前进角为50°～60°，后退角0°～10°。

式(6)、式(7)就构成了孔隙网络准静态两相流动模拟的主要微观填充机理，初次排驱较简单，但对于二次吮吸过程，多个微观机理相互之间是竞争关系，需要按照公式中所计算的各类填充临界毛细管力的大小进行排序，较大的临界毛管力对应的填充机理优先发生。但由于孔隙喉道输入巨大，每一步填充一个孔或喉是不现实的，为获取更高的效率，每一个毛细管力变化内，有很多孔和喉均有填充发生，但填充发生的早晚影响着协同式孔隙填充发生过程中n的值，为此提出：首先，依据式(6)确定所有卡断填充发生的喉道并令其填充，继而，依据式(7)判断的协同式填充是否发生，但对于n=0时，也需要额外按照n=1判断孔隙内油是否为死油，最后，与已被水填充的孔隙的喉道也令其均被水填充。

此外，一个重要的问题是，油相作为非润湿相会以油团的形式被圈闭在孔隙网络中，形成残余油，油团可能是一个孔[式(7)，n为0时]，也可能是若干个孔隙喉道的连通体，此过程需要高效算法实时监测不可动油团的形成[13]。通过采用先进的图论算法，借助NetworkX[19]高效实现了对孔隙网络连通关系的动态监测。其主要思想为：①增加 Outlet 节点与出口处所有孔节点连接；②找到与Outlet连接的所有含油的孔节点，即为可动油；③所有含油孔节点去掉可动油孔节点即为不可动的死油团。

1.3 油湿及混合润湿油藏水驱油过程模拟

当储层为油湿时，水驱油过程中，毛细管力为阻力，此过程类似于初次排驱过程，区别在于存在水膜，主要的微观充填机理有带油膜活塞式充填和不带油膜活塞式充填。

(8)

(9)

此外，对于混合润湿岩心水驱油过程的模拟。认为孔隙半径较大的孔隙为油湿孔，与之连接的喉道也为油湿，给定油湿孔隙体积所占比，标记油湿孔隙和喉道。对于水湿的孔隙和喉道所组成的孔隙网络按照2.2节开展模拟计算，此时保持油湿孔隙网络均为存在水膜的油相充填，并按照整体的孔隙网络计算相对渗透率曲线和毛细管力，继而按照上述算法对油湿的孔隙网络开展模拟计算，最终获得整体的相对渗透率曲线和毛细管力曲线。

2 宏观物性参数计算方法

首先开展单相流的模拟计算，获取岩心的绝对渗透率。具体地，计算每个孔隙和喉道的传导系数T，计算公式为

(10)

式(10)中：T为传导系数，m3；L为岩芯长度，m。

对于孔隙i与孔隙j之间的传导系数可按照调和平均获得，可表示为

(11)

式(11)中：Ti为孔隙i的传导系数，m3；Tj为孔隙j的传导系数，m3；Tij为喉道的传导系数，m3。

孔隙尺度模拟均假设流体不可压缩，因此对于任意一个孔隙与周围所连接孔隙之间的流入和流出量之和为零，以此可建立以孔隙压力为未知量的线性方程组。基于图论的处理方式，利用孔隙喉道之间的连通关系，建立邻接矩阵，邻接矩阵的系数即为传导系数。添加邻接矩阵的对角线为每一行矩阵累加和的负数，并且根据边界孔隙编号，将边界孔隙对应的矩阵进行处理，设定入口和出口压力值。将所建立的矩阵变为稀疏矩阵，利用大规模稀疏矩阵的求解算法，获得每一个孔隙的压力。继而根据压力场和传导系数，获取出口端面的流量q，利用达西定律获得绝对渗透率，即

(12)

式(12)中：q为流量，m3/s；μ为液体黏度，Pa·s；pin、pout分别为入口和出口处的压力，Pa。

对于两相流过程，对于每个时间步，根据饱和度和截面内油水的分布形式，确定计算每个孔隙和喉道的油相和水相的传导系数[20]

(13)

To=T(1-Sw)2

(14)

式(13)中：Tw为水相的传导系数，m3；To为油相的传导系数，m3；C(G,θ)为润湿相额外附加流动阻力的修正系数，具体值参考文献[21]。

最后按照与单相流计算相同的方式，计算油相和水相的压力场，获取油相渗透率ko、水相渗透率kw，即可获得油、水相对渗透率

而对应的整体含水饱和度可通过孔和喉的体积加权获取。

3 实例应用

3.1 Berea砂岩岩心

kr为相对渗透率图2 油驱水(初次排驱)过程相对渗透率曲线Fig.2 Relative permeabilities of oil injection (primary drainage) process

读入已获得的经典Berea砂岩的孔隙网络StateOil格式数据[9]，按照上述计算方法，对于油驱水(成藏)、水湿岩心水驱油以及油湿岩心水驱油过程进行了模拟。本代码效率较高，普通笔记本电脑模拟计算时间少于5 min，得到图2、图3所示的模拟计算结果，与室内实验测量结果[22]吻合率高，证实了所建立孔隙网络流动模拟方法的有效性和准确性。由图3可知，水湿岩心相对渗透率曲线残余油饱和度高，且水相相对渗透率端点值小于0.2；而油湿岩心相对渗透率曲线残余油饱和度低，且对应的水相端点值很高，接近0.8。这些相对渗透率曲线特征是与前人的研究是一致的[23]。

图3 岩心水驱油过程的相对渗透率曲线Fig.3 Relative permeabilities of water flooding process

在此研究基础之上，对接开源软件Paraview实现了两相流动的三维可视化表征。图4为水湿和油湿情况下水驱剩余油分布，其中孔隙以球体的形式表征，而喉道以棍棒的形式表征，其颜色表征着流体的饱和度，红色为油相，蓝色为水相。图4可证实水湿岩心剩余油主要以孤立油团的形式存在，且主要分布于孔喉比较大的孔隙内，而油湿岩心水驱剩余油主要分布于小孔喉的区域，而大孔动用程度高，与水湿情况刚好相反。

红色为油相；蓝色为水相图4 水驱后剩余油分布的三维可视化表征Fig.4 Three dimensional visualization of water-flooded rock

3.2 吉木萨尔页岩

在Berea砂岩两相流动模拟结果与实验验证的基础上，利用FIB-SEM实验实现了对吉木萨尔页岩的三维孔隙空间数字化表征。其分辨率为20 nm×20 nm×20 nm，数据空间大小为500×500×500。应用快速傅立叶变换算法，去除了三维图像中的噪声伪影，继而通过设置阈值灰度值对孔隙空间进行图像分割，并采用最大球体算法提取页岩的三维孔隙网络。具体的页岩三维孔隙空间和孔隙网络的Paraview可视化如图5所示。

图5 吉木萨尔页岩孔隙空间的三维可视化及孔隙网络模型可视化表征Fig.5 Three dimensional visualization of pore space and pore network of Jimsar shale

吉木萨尔页岩孔隙网络共有13 419个孔隙和31 393个喉道，连通孔隙度为16.35%，与实验室内气测孔隙度19%接近。孔隙率的损失是由于小于20 nm的孔和喉道无法识别造成的。孔隙和喉道的平均半径分别为29.75 nm和19.13 nm。尽管有几个相当大的孔和喉，但整体平均孔喉比为2.29。应用所提出的孔隙网络方法，计算出x、y和z方向的达西渗透率为0.004～0.006 mD，这也是页岩油藏在地层条件下的合理值。页岩储层通常认为是混合润湿，此处采用1/2水湿和1/2油湿的混合润湿条件，将较大/较小部分的孔隙和喉道设置为油湿，其余为水湿情况，模拟计算得到的油水相对渗透率曲线如图6所示。

图6 吉木萨尔页岩混合润湿条件下油水两相相对渗透率曲线Fig.6 Oil-Water relative permeability curves of mixed-wetted Jimsar shale

由图6可知，两种混合润湿性情况均表明残余油饱和度较高(即接近40%)，且两相流区域狭窄，这意味着页岩油的采收率远低于常规高渗透油藏。对于较大孔隙油湿情况，水相相对渗透率极低，且端点值小于0.01，这意味着页岩基质中的水几乎不会流动，这与致密储层“渗透率圈闭”现象是一致的[24-26]。对于较小孔隙油湿情况，随着含水饱和度的增加，水相相对渗透率增加得更快，端点值在0.3左右。两者的差异表明，较大的孔隙具有主要的流动能力，当较大的孔隙为油湿时，水首先充满水湿的较小孔隙，然后充满油的较大孔隙被卡断形成死有团。然而，当较大的孔隙为水湿时，水首先填充孔隙，然后延伸至较小的油湿孔隙。油湿孔隙的分布对油水流动机制起着重要作用，未来的研究需要建立更大尺度的孔隙网络模型和更精细的复杂润湿性表征，使其更能代表页岩地层，这是未来工作的研究方向。

4 结论

基于实际数字岩心抽提的孔隙网络模型，利用四角星形截面流体微观赋存机理，结合高性能图论算法和开源可视化表征软件，考虑了活塞式驱替、卡断填充、协同式填充等孔喉微观填充机理，建立了孔隙网络两相流动模拟方法。实现了对孔隙网络内非润湿相连通关系的动态监测，模拟了初次排驱和水驱油过程，基于经典Berea砂岩的孔隙网络，计算得到了不同润湿性条件下油水的相对渗透率曲线，对比实验结果验证了模拟结果的有效性和准确性，可视化表征结果方便的表征了不同润湿条件下剩余油的分布规律。得出如下主要结论。

(1)证实了水湿岩心剩余油主要以孤立油团的形式存在，且主要分布于孔喉比较大的孔隙内，而油湿岩心水驱剩余油主要分布在小孔喉的区域，而大孔动用程度高，与水湿情况刚好相反。

(2)对吉木萨尔页岩储层油水两相流动进行了初步探讨，结果表明页岩两相流区域狭窄，残余油饱和度高，水相存在“渗透率圈闭”现象，而且明确了油湿孔隙的分布对油水流动机制起的重要作用。