CN113468829B - 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法 - Google Patents

Abstract

本发明涉及油气田开发技术领域，涉及一种基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，包括：1)由核磁共振实验与恒速压汞实验获取锥状孔喉半径分布；2)获取锥状孔喉长度和配位数；3)建立孔隙网络模型；4)在孔隙网络模型中引入非稳态渗流理论，同时根据孔隙内流体流动、界面移动和孔隙压力扩散过程，研究孔隙型介质内部水驱油两相非稳态渗流，形成非牛顿两相流体驱替模拟方法，将动态网络模拟算法与非稳态渗流理论相结合模拟非牛顿两相流体驱替过程。本发明能准确反映非牛顿两相流体在多孔介质孔喉内部流动过程中，流体与孔喉压缩系数、流体黏度等随压力传播而变化的非稳态渗流特征。

Description

基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法

技术领域

本发明涉及油气田开发技术领域，具体地说，涉及一种基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法。

背景技术

地下原油储量大，用途广，是现代工业的“血液”，更是国家生存和发展不可或缺的战略资源，如何最大程度的将地下能源采出来转变成经济发展的动力显得尤为重要，其中的根本就是如何提高油田采收率。经过多年来的研究，常规油藏的开发方式已经趋于成熟，研究人员大都采用符合储层岩石孔隙微观结构特征的孔隙网络模型进行油藏数值模拟研究，且该类技术已经日趋成熟，但在孔隙网络模型中考虑开采过程中油藏非稳态渗流机理研究还比较少。

目前，国内外许多研究人员通常采用连续介质理论研究多孔介质的多相渗流，但由于黏滞力、毛管力在孔隙尺度上具有非连续性，且当流体介质为非牛顿流体时的渗流规律如何变化尚不明确；实际流体在储层内部通常具有可压缩性，这与流体能瞬时从入口传到出口的稳态渗流理论相悖。在油藏中，参与渗流的很多流体都是非牛顿流体，如提高原油采收率工艺中所用到的微乳液、聚合物溶液、各种酸化液和压裂液均为非牛顿流体。相比于牛顿流体，非牛顿流体大多由高分子物质所组成，而高分子物质往往会在分子力的作用下形成各种盘绕的网状或链状结构，引起流动过程中黏度变化，因此，牛顿流体与非牛顿流体在宏观上的明显差异为在不同的剪切作用下黏度不同，故而牛顿流体与非牛顿流体的流变特征(流体受到外力作用时发生流动和变形的性质叫做流变性)存在很大的区别，在多孔介质内部流动所表现出来的渗流性质也有很大差异，而且研究也比较困难。

现国内外基于格子Boltzmann方法的非牛顿流体流动研究较多，但基于孔隙网络模型的两相非牛顿流体非混相驱替模拟研究仍较少，且大都无法准确反映两相流体在多孔介质孔喉内部流动过程中，无法描述流体与孔喉压缩系数等随压力传播而变化的非稳态渗流特征，该类数值模拟结果往往并不十分可靠，不能为现场生产提供科学指导。因此，提出一种适用于非牛顿流体的油藏孔隙网络模型非稳态渗流模拟方法是十分必要的。

发明内容

本发明的内容是提供一种基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其能够克服现有技术的某种或某些缺陷。

根据本发明的一种基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其包括以下步骤：

(1)、由核磁共振实验与恒速压汞实验获取锥状孔喉半径分布；

(2)、通过微CT扫描得到锥状孔喉长度和配位数；

(3)、根据锥状孔喉半径分布、孔喉长度和配位数，建立孔隙网络模型；

(4)、在孔隙网络模型中引入非稳态渗流理论，同时根据孔隙内流体流动、界面移动和孔隙压力扩散过程，研究孔隙型介质内部水驱油两相非稳态渗流，引用非牛顿流体，形成非牛顿两相流体驱替模拟方法，将动态网络模拟算法与非稳态渗流理论相结合模拟非牛顿两相流体驱替过程。

作为优选，步骤(1)中，建立SC模型，并在SC模型的基础上构建无序孔隙网络模型；通过核磁共振实验得到岩心孔径分布T2谱，结合与恒速压汞实验，并根据预设定量关系将所述T2谱转换得到孔喉半径分布频率，根据所述孔喉半径分布频率拟合得孔喉分布函数。

作为优选，步骤(2)中，微CT扫描得到配位数与喉道长度包括以下步骤：

(2.1)利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品和射线源及探测器进行360°的相对旋转，采集上千帧角度的数据，然后利用计算机断层扫描成像重构方法进行3D重构，从而得到样品内外结构的高分辨3D数据及影像，从而完成利用微CT扫描实验获取CT实验数据；

(2.2)进行CT数据分析，通过图像不同灰度进行物质区分，灰度值低的区域代表物质密度低，参考灰度曲线中孔隙的灰度值进行阈值划分，在图像中将孔隙进行单独的提取分离；

(2.3)利用最大球算法，区分数字岩心三维图像中的孔隙、喉道所占空间及连通性，提取相应的孔隙、喉道结构网络模型，同时运用数理统计方法可实现对孔喉长度、孔喉体积、孔喉比、配位数、形状因子的定量提取，得到研究岩石孔喉表征的参数；

其中，喉道长度l可由下式计算：

l＝D-R₁-R₂；

式中，R₁，R₂分别为该喉道所连接两个孔隙的半径；D为两孔隙中心点的实际坐标距离。

作为优选，步骤(3)中，建立孔隙网络模型的具体步骤如下：

(3.1)确定网络模型的大小与节点数，构建一个X×Y×Z的三维简单立方体网格；每个节点代表一个孔隙，节点与节点之间由喉道相连；由此建立起来的网络中每个代表孔隙的节点周围都有六个喉道相连，每个喉道两端与两个孔隙相连；其x方向、y方向和z方向各节点之间的间隔距离为喉道长度l，各方向节点数为nx,ny,nz，模型的边长为L_x＝(nx-1)×l、L_y＝(ny-1)×l、L_z＝(nz-1)×l；

(3.2)计算出网络模型中各节点的坐标；计算式为：

(x，y，z)＝[(i-1)l，(j-1)l，(k-1)l]；

式中i、j、k分别为x方向、y方向和z方向的节点序号；

(3.3)通过中心点位移中心节点坐标生成随机网络；

按照如下式移动各节点坐标(x，y，z)，实现各节点坐标在半径为0.5l的球面区域内随机移动，rand()为随机函数；

(3.4)在程序中设定概率为p的概率函数；

p＝z/z_max，z为配位数，表示孔喉间的连接概率，z_max为SC模型具有的最大配位数，其值为6；

如当连通概率p为66.7％时，此时对应配位数为4，在rand()函数生成的1～100整数中有66.7％的概率小于50，有另外33.3％的概率大于50；因此，将小于50的数值的管束与节点连接，并将其索引值赋为1，不满足条件的赋值为0，从而为每个管束赋予半径属性；

(3.5)为模型赋予孔喉半径属性与喉道长度属性，首先利用C++语言，借助科学计算库Eigen，生成2个与孔喉数量相等长度的单位向量，取名为半径向量与喉道长度向量，此时半径向量与喉道长度向量中的所有元素值均为1，表示模型中管束的半径值与喉道长度均为单位1，再通过获取的孔喉分布函数为半径向量中的每个元素赋值，这样就产生了符合孔喉半径大小的“半径向量”，完成孔喉半径赋值，从而生成一个符合岩心孔喉半径频率分布的三维张量数据体；通过对喉道长度向量中每个元素赋喉道长度值；

若孔隙网络模型孔喉形状为圆管状，则圆管状泊肃叶定律表达形式为：

孔隙网络模型孔喉形状为锥状时，则锥状孔喉的泊肃叶公式推导过程如下：

设流体通过管道微元体，通过长度为dx时的流体流动状态是层流的，并且遵循雷诺方程，此时泊肃叶定律表示为：

其中p(x)是在x，且平滑函数r(x)＝r_t+x(r_p-r_t)/l_c处的压力，l_c为锥状管道长度，r_p为孔隙半径，r_t为微元dx所在部分的喉道半径，r_pt为孔喉比，因此：

因此可以得到锥状孔喉泊肃叶公式为：

可写为：

其中：

当孔喉为均匀的管柱，孔喉比为1，α＝1，此时的锥形管束泊肃叶公式等于圆形管柱泊肃叶公式。

作为优选，步骤(4)中，非稳态非牛顿两相流体模拟方法为：

孔隙网络模型中，管束内两相流体共存时非润湿流体与润湿流体界面引起的毛细管压力p_cij通过Young-Laplace方程计算：

p_cij＝-2γcosθ/r_ij；

其中，r_ij为管束两端节点i和节点j之间管束中两相界面处的半径，γ为两相界面张力，θ为润湿角；

毛管压力为弯液面在管束中的位置的函数，称为动态毛管压力p_cij，因此，将Young-Laplace改为：

其中X_ij是与凹液面位置有关的无量纲数，即凹液面所在位置横坐标除以整个孔隙网络的长度，当管道中只有单相流体时，p_cij＝0；且管束内两相共存时，当弯液面处于管束末端时，p_cij＝0，在管束中间时具有最大值p_cij＝4γ/r_ij；

孔隙网络模型中；当存在两相弯液面(即两相共存)时，管内流量计算公式为：

其中其中r_ij为管束半径，p_i和p_j分别为管束两端节点i和节点j的孔隙压力，μ_eff为管束内有效黏度，其计算方法为μ_eff＝μ_IX_ij+μ_D(1-X_ij)，其中μ_I为非润湿相(驱替相)流体粘度，μ_D为润湿相(被驱替相)流体黏度，初始孔隙网络模型中全部充填油相，作为被驱替相采出，模型润湿角取180，水相为非润湿相作为驱替相注入，当管束内只有非润湿相流体时，μ_eff＝μ_I，p_cij＝0，当管束内只有润湿相流体时，μ_eff＝μ_D，p_cij＝0；

对非牛顿流体中的宾汉模型与Ellis模型进行研究：

(1)宾汉流体；当剪切应力超过某一静态剪切应力时，流体发生流动，其本构方程为：

τ＝τ₀+μ_bγ；

其中，τ₀为屈服应力，μ_b为宾汉流体黏度；

当管束内只有宾汉流体时，流量计算公式为：

其中r_ij0为岩心物理模型半径，在模型中，假设模型中所有管束的r_ij/r_ij0均相同，则对宾汉流体来说，上式可写为

其中β表示与非牛顿流体流变性有关的参数，在模拟过程中，被驱替相只改变了黏度，每次迭代计算过程中，随压力变化时，黏度μ_b通过计算公式μ_bnew＝(μ_b/1-β)更新。

当管束内水与宾汉流体两相共存时，μ_eff的计算中，μ_D＝μ_b；若β＝0，则被驱替相为牛顿流体；模拟过程中，更改β，即可更改宾汉流体性质；

(2)Ellis模型；每次迭代计算过程中，随压力变化时，Ellis流体的黏度μ_e通过如下计算公式更新：

其中，μ_o是低剪切应力时的黏度，τ_1/2是μ_o＝μ_o/2时的剪切应力，α是经验常数；

当管束内只有宾汉流体时，流量计算公式为：

当管束内水与Ellis流体共存时时，流量计算公式为：

此时μ_D＝μ_e；模拟过程中，更改τ、α和μ_o，即可更改Ellis流体性质；

根据质量守恒定律构建类似AP＝B的矩阵方程，通过逐次超松弛迭代法求解，同时需要根据时间步长Δt校正；驱替过程一直进行直到整个孔隙网络空间被侵入流体占满；然后重新计算并更新压力场和饱和度，进行下一步计算；整个过程保持注入流量Q不变。

本发明提出了一种适用于非牛顿流体的油藏孔隙网络模型非稳态渗流模拟方法，能准确反映非牛顿两相流体在多孔介质孔喉内部流动过程中，流体与孔喉压缩系数等随压力传播而变化的非稳态渗流特征，能为现场生产提供科学指导。

附图说明

图1为实施例1中一种基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法的流程图；

图2为实施例1中孔喉形状抽象为锥状的管束模型；

图3为实施例1中孔隙空间中只有活塞式驱替的示意图；

图4为实施例1中由压汞与核磁实验获取的孔喉分布图；

图5为实施例2中SC网络模型图；

图6为实施例2中水驱宾汉流体流量图；

图7为实施例2中水驱Ellis流体流量图。

具体实施方式

为进一步了解本发明的内容，结合附图和实施例对本发明作详细描述。应当理解的是，实施例仅仅是对本发明进行解释而并非限定。

实施例1

如图1所示，本实施例提供了一种基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其技术方案如下：

1.由核磁共振实验与恒速压汞实验获取孔喉半径分布

1.1核磁共振实验

核磁共振实验原理可以简述为原子核的自旋能级在外加磁场环境中发生相应的分裂并产生共振现象。原子核外的电荷围绕原子核心不停旋转，就产生了一个具有强度和方向的矢量磁场，如果没有外加磁场，单个核磁矩的方向是随机的，则在宏观上不表现出磁性。当检测对象置于外加磁场中时，核磁矩将向外加磁场缠绕的方向进动，若再增加一个外加交变电磁场，它将从低能量状态吸收能量并跳到高能量状态。一旦进入高能状态后取消外加射频，原子核就会逐渐从高能态向低能态恢复，这个过程就是弛豫过程，这一过程所需要的时间就叫弛豫时间。弛豫过程分为纵向弛豫过程和横向弛豫过程，对应纵向弛豫时间T₁和横向弛豫时间T₂，弛豫过程又由体积弛豫(b)、表面弛豫(s)、和扩散弛豫(d)三部分构成。所以，对于岩石中油相和水相的纵向弛豫时间T₁和横向弛豫时间T₂则可以表示为：

核磁共振的观测对象是氢核，在评价储层岩石时，岩石骨架中的氢核弛豫时间极短，且横向弛豫过程比纵向弛豫过程衰减更快，在核磁共振仪器采集到回波信号之前，其横向弛豫过程基本已经衰减完毕完毕，所以，将横向弛豫时间T₂作为主要的评价参数，可以基本排除岩石骨架的信号。对于饱和流体的岩心，水相主要表现出表面弛豫，油相主要表现出自由弛豫及扩散弛豫，岩心饱和流体后核磁共振扫描捕捉到多种弛豫信号的叠加，最终核磁共振观测到的回波串信号可用多指数衰减规律表示，如式3：

式中，M(t)为t时刻观测到的回波信号幅度；M_i为第i种弛豫在零时刻的回波信号幅度；T_2i为第i种弛豫的横向弛豫时间，ms。

在使用核磁共振技术进行孔喉结构表征时，常以完全饱和水的岩样横向弛豫时间来表征岩样的完整孔喉分布，此时横向弛豫时间主要由表面弛豫决定，如式4：

其中，ρ₂为表面弛豫率；V为孔隙体积，μm³；S为孔隙表面积，μm²。

式中S/V还可以表示为孔隙无因次形状因子(F_s)与孔隙半径(r)的比值。球形孔隙F_s＝3，圆柱形孔隙F_s＝2，则孔喉半径r与弛豫时间T₂的关系如下式，所以岩石孔隙中流体的横向弛豫时间与孔隙大小、形状相关，孔隙越小，对应T₂值越小。

对某一给定的岩样，ρ₂和F_s近似为常数。因此，可使用一个恒定的转化系数C表示1/ρ₂F_s，即：

由式6可知，T₂弛豫时间与水赋存孔隙的孔径大小存在对应关系，即r＝CT₂，如果计算出转化系数C，则可完成横向弛豫时间转换为孔隙半径，获得孔喉半径分布。

1.2恒速压汞实验

恒压压汞是一种常规的压汞测试方式，其原理是在某一恒定压力条件下将非润湿相汞注入岩心孔隙并计量，依次升高注汞压力注汞并多次计量，最终获得进汞压力和饱和度之间的关系，即毛管压力曲线。当注汞压力达到最大值后，依次减小注汞压力直至无汞退出，则得到退汞压力曲线，用于分析岩石表面润湿性及孔隙介质中驱油机理。根据注汞压力计算不同压力条件下对应的孔喉半径和孔喉大小，得到岩样孔喉大小分布，较低注汞压力对应测量较大孔隙，当注汞压力依次升高，汞开始进入更小的孔喉，当最小的孔喉被汞充填时，达到汞饱和度的最大值。根据毛管压力理论公式(式7)，由注汞压力可直接计算出对应的孔隙大小。

式中，p_c为毛细管压力，单位为MPa；σ为界面张力，单位为mN/m²；θ为静态接触角，单位为°；r为毛细管半径，单位为μm。

目前，计算转化系数的方法主要为基于核磁共振T₂谱构造伪毛管压力曲线，或基于毛管压力曲线构造伪T₂谱。这两种方法在本质上是一致的，即转换为同一坐标条件来对比刻度，从而求得转化系数，以获取岩样完整的孔隙分布信息。本实施例利用恒速压汞实验获取平均孔喉半径，并通过与核磁共振实验平均T₂谱中值之比获取转换系数c。且所述转化关系为r_m＝cT_2m，r_m为第m个孔喉的半径，T_2m为T₂谱的第m个幅度值，m为正整数；最后根据孔喉半径分布频率拟合得到孔喉分布函数，具体包括：通过随机对数正态分布随机函数拟合，其中，通过拟合所述孔喉半径频率分布得到数学期望μ和标准偏差σ，x表示孔喉半径。

2.微CT扫描实验统计配位数与喉道长度

利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品和射线源及探测器进行360°的相对旋转，采集上千帧角度的数据，然后利用计算机断层扫描成像重构方法进行3D重构，从而得到样品内外结构的高分辨3D数据及影像，从而完成利用微CT扫描实验获取CT实验数据(主要是图像数据)。利用计算机断层扫描成像重构方法进行3D重构，从而得到样品内外结构的高分辨3D数据及影像，CT数据分析是通过图像不同灰度进行物质区分，灰度值低的区域代表物质密度低，参考灰度曲线中孔隙的灰度值进行阈值划分，在图像中将孔隙进行单独的提取分离。利用最大球算法，区分数字岩心三维图像中的孔隙、喉道所占空间及连通性，提取相应的孔隙、喉道结构网络模型，同时运用数理统计方法可实现对孔喉长度、孔喉体积、孔喉比、配位数z、形状因子等孔隙结构的定量提取，得到研究岩石孔喉表征的参数。

其中，喉道长度l可由式8计算：

l＝D-R₁-R₂ (8)

式中，R₁，R₂分别为该喉道所连接两个孔隙的半径，单位为μm；D为两孔隙中心点的实际坐标距离，单位为μm。

3.建立孔隙网络模型

基于孔隙网络模型的数值模拟方法一般假设岩石内部孔喉具有一定的形状，且可以通过岩心分析实验准确获取这类资料，模型复用性高，相比实验，可变条件(流速、压力等)宽泛，且稳定的孔隙网络模拟方法可无限重复使用，具有强大应用价值与经济效益。以SC网络模型为基础通过中心点位移生成无序结构网络，利用C++语言建立盘状网络模型的过程，其具体步骤如下：

(1)确定网络模型的大小与节点数(或网格数)，构建一个X×Y×Z的三维简单立方体网格。每个节点代表一个孔隙，节点与节点之间由喉道相连。由此建立起来的网络中每个代表孔隙的节点周围都有六个喉道相连，每个喉道两端与两个孔隙相连。其各方向(即x方向、y方向和z方向)各节点之间的间隔距离设置为l(l为喉道长度)，各方向节点数设置为nx,ny,nz，模型的边长为L_x＝(nx-1)×l、L_y＝(ny-1)×l、L_z＝(nz-1)×l。

(2)计算出网络模型中各节点的坐标。计算式为：(x，y，z)＝[(i-1)l，(j-1)l，(k-1)l]，式中i、j、k分别为x方向、y方向和z方向的节点序号，取值分别为1，2，3…。

(3)通过中心点位移中心节点坐标生成随机网络。

按照如公式(9)移动各节点坐标(x，y，z)，实现各节点坐标在半径为0.5l的球面区域内随机移动，rand()为随机函数。

(4)在程序中设定概率为p的概率函数(p＝z/z_max，z为配位数，表示孔喉间的连接概率，z_max为SC模型具有的最大配位数，其值为6))，如当连通概率p为66.7％时(此时对应配位数为4)，在rand()函数生成的1～100整数中(每个整数对应一个喉道)有66.7％的概率小于50，有另外33.3％的概率大于50。因此，将小于50的数值的管束与节点连接，并将其索引值赋为1(索引值为1的管束，按照孔喉分布函数来给各管束半径赋值)，不满足条件的赋值为0(半径值设为0，不分配半径)，从而为每个管束赋予半径属性。需要注意的是，三维网络模型中存在一个临界配位数z_c，这是模型能够连通的最小配位数，一般情况下三维模型z_c为1.5。

(5)为模型赋予孔喉半径属性与喉道长度属性，首先利用C++语言，借助科学计算库Eigen，生成2个与孔喉数量相等长度的单位向量(取名为半径向量与喉道长度向量)，此时半径向量与喉道长度向量中的所有元素值均为1，表示模型中管束的半径值与喉道长度均为单位1，再通过获取的孔喉分布函数为半径向量中的每个元素赋值，这样就产生了符合孔喉半径大小的“半径向量”，完成孔喉半径赋值，从而生成一个符合岩心孔喉半径频率分布的三维张量数据体(称为基础孔隙数据体，通过该数据体即可绘制三维孔隙介质模型)；通过对喉道长度向量中每个元素赋喉道长度值(本实施例认为所有喉道长度相同，因此该向量为常数向量)。

上述孔隙网络模型假设孔喉形状为圆管状，且圆管状泊肃叶定律通常表达形式为：

考虑到真实岩心孔喉的复杂性，并且综合考虑模拟运算效率与孔隙网络模型精度后，本次采用网络模拟采用锥状孔喉构建网络模型，锥状孔喉配置结构示意图如图2，适用于锥状孔喉的泊肃叶公式推导过程如下：

设流体通过图2中的管道微元体，通过长度为dx时的流体流动状态是层流的，并且遵循雷诺方程，此时泊肃叶定律表示为：

其中p(x)是在x，且r(x)＝r_t+x(r_p-r_t)/l_c(又称平滑函数)处的压力，l_c为锥状管道长度，r_p为孔隙半径，r_t为微元dx所在部分的喉道半径，r_pt为孔喉比，因此：

因此可以得到锥状孔喉泊肃叶公式为：

对比式(10)，式(13)可写为：

其中：

式(14)中，当孔喉为均匀的管柱，孔喉比为1，α＝1，此时的锥形管束泊肃叶公式等于圆形管柱泊肃叶公式。相比圆柱状孔喉结构，锥状管束的优点在于：(1)锥状孔喉泊肃叶公式可以通过转换系数α从圆管泊肃叶定律计算得到，计算更方便；(2)可以将真实岩心孔喉比r_pt值直接引入孔隙网络模拟计算中，相比于孔隙-喉道-孔隙串联模型，该方法更简单直接且具有可靠性；(3)模拟过程中通过转换系数α可以在保证模型情况更符合真实情况的同时大大提升计算效率。

4.动态网络模拟方法

4.1非稳态单向流模拟方法

在真实岩心中，由于流体的粘性作用，流体质点粘附在物体表面上，形成流体不滑移现象(即相对速度为零)，因而产生摩擦阻力和能量耗散。因此假设孔隙网络中的流体流动遵循能量耗散最低原则，且粘性流动仅在流体窜流分支方向发生，流动过程中孔隙网络模型遵循的质量守恒定律通过基尔霍夫定律描述，即流入的流体体积等于流出的流体体积，由此将真实岩心基质流动简化为孔隙网络模型流动。

根据基尔霍夫定律，引入无流动边界条件，解得各节点的流动压力，由此求得各截面的平均流速。模拟过程中整体流动方向设定为水平方向，流体为单相，模型中对于单独的节点i和j，设两点压力分别为p_i和p_j，两点间连接孔道的半径和长度分别为r_ij和l_ij，流体粘度为μ，水力传导系数为则两节点间的体积流量用泊肃叶原理可以表达为：

传统的孔隙网络模拟中，根据基尔霍夫定律，两节点间的总流量为零。则流体在孔隙网络中流动满足拉普拉斯方程：

然而，实际渗流过程中，通常认为孔隙内部流体和岩石骨架具有微可压缩性。若考虑流体和岩石微可压缩性，则孔隙网络模型中的单相非稳态渗流方程为：

其中，其中为Hamilton算子；g表示孔隙网络模型中具有喉道连通的节点间的流体传导率，p表示孔隙网络模型的节点压力，C_t表示孔隙网络模型的节点的综合压缩系数，表示节点压力对时间的偏导数。

孔隙网络模型中的最小单元(微元体)为节点和管束，对应的流体质量(体积)流量满足泊肃叶定律。因此，通过泰勒展开以及隐式有限差分法，将方程(12)应用到孔隙网络模型中，并考虑每个节点中注入(记为正)或流出(记为负)的流量(q_IOi)后得到：

其中，g_ij表示节点i和j之间的管束的传导率，p_i和p_j分别表示节点i和与之相邻的节点j的节点压力，V_bi表示节点i和j之间的管束的体积，dp_i/dt表示不同时间步长时压力的变化(时间不长确定方法见2.3节)，且i为小于等于N的正整数，N为孔隙网络模型中全部节点数，c为小于等于节点i的配位数的正整数。

根据上式，遍历孔隙网络中的所有节点可构建出单相流线性方程组或矩阵方程：

通过转化得到矩阵方程AP＝B，其中A为对称、正定、主对角占有的稀疏方阵，且其中任意元素A_i可表示为：1)当A_i为非对角线元素时，A_i＝g_ij，2)当A_i为非对角线元素时，上标n表示当前时刻；P为由各节点压力构成的向量[P]＝[p₁,p₂,p₃,···p_N]^T；B为向量[B]＝[B₁,B₂,B₃,···B_N]^T，且向量B中任意元素B_i可表示为：/>

可使用超松弛迭代法(SOR)求解矩阵方程AP＝B，迭代公式如下:

式中下标ij表示连接第i节点的第j个节点，α是松弛因子(经试算取1.75)。在得到压力场后，即可算得整个孔隙网络的流量通量q(m³/s)。

4.2非稳态非牛顿两相流体模拟方法

孔隙网络模型的多相流数值模拟方法有：(准)静态网络模拟和动态网络模拟。(准)静态网络模拟只考虑毛管力对渗流的影响，忽略流体的黏性压降和时间等影响因素。动态网络模拟技术可研究孔隙尺度多相流的瞬变现象，可同时考虑毛管压力、黏滞力、重力和时间等对多相流的影响。

孔隙网络中的多相流流动模拟存在以下假设条件：(1)所有的流体都被认为包含在孔隙管道和节点中，但是所有的压降都发生在节点之间的管道中；(2)孔隙网络中两相流体间只存在一个界面；(3)网络中有两种不可混溶的流体流动；(4)两种流体在管口界面上的毛细管压差与管半径成反比；(5)流体流量可以用泊肃叶方程计算；(6)孔隙空间中只有活塞式驱替(图3)。

孔隙网络模型中，管束内非润湿流体与润湿流体界面引起的毛细管压力p_cij通过Young-Laplace方程计算：

p_cij＝-2γcosθ/r_ij (21)

其中，r_ij为两相界面(图3中弯液面)处半径，γ为两相界面张力，θ为润湿角(其范围为0-π/2)。

在本发明中，我们假设管束为沙漏形状，弯液面处的有效半径r_i遵循平滑函数。因此，毛管压力也为弯液面在管束中的位置的函数(称为动态毛管压力p_cij)，因此将Young-Laplace改为：

其中X_ij是与凹液面位置有关的无量纲数，其值定义为管束内界面位置与管束长度之比(0≤X_ij≤1)，当管道中只有单相流体时，p_cij＝0；且管束内两相共存时，当弯液面处于管束末端时，p_cij＝0，在管束中间时具有最大值p_cij＝4γ/r_ij.

孔隙网络模型中，管束内为单相流体时，管束内流体采用公式(15)计算，其中粘度为该相流体黏度；当存在两相弯液面时(即管束内两相共存时)，管内流量计算公式为：

其中μ_eff为管束内有效黏度，其计算方法为μ_eff＝μ_IX_ij+μ_D(1-X_ij)，其中μ_I为非润湿性流体粘度，μ_D为润湿相流体黏度，由于本发明为水驱油，因此初始孔隙网络模型中全部充填油相，作为被驱替相采出，模型润湿角取180(完全油湿)，水相为非润湿相作为驱替相注入。

实际上，许多种类型的油均为非牛顿流体(如稠油，含矿物微粒原油等)，这类油品在温度降低时，其性质逐渐偏离牛顿流体(如粘度变化等)。因此，多孔介质中的大量油品是非牛顿流体。本发明假设每个喉道内的驱替流体都是牛顿流体，且被驱替的流体是非牛顿流体，本发明主要对非牛顿流体中的宾汉模型与Ellis模型进行了研究。

(1)宾汉流体(也称塑性流体)。当剪切应力超过某一静态剪切应力时，流体发生流动，其本构方程为：

其中，τ₀为屈服应力，μ_b为黏度，为应变率张量。

当管束内只存在宾汉流体时，流量q_ij计算公式为：

其中r_ij0为岩心物理模型半径，在模型中，假设模型中所有管束的r_ij/r_ij0均相同，则对与宾汉流体，式(24)可近似写为

其中β表示与宾汉流体流变性有关的参数，式(42)表示，在模拟过程中，随时间步长变化，被驱替相(非牛顿流体)黏度发生改变，改变后的黏度计算公式为n+1为下一时刻，n为当前时刻。

当管束中存在两相流体时，此时，

且μ_eff的计算中，μ_eff＝B_wμ_wX_ij+B_o(μ_o/(1-β))(1-X_ij)，μ_o为当前时刻的流体黏度，随时间步变化而根据公式为变化。若β＝0，则被驱替相为牛顿流体。模拟过程中，更改β，即可更改宾汉流体性质。

(2)Ellis模型。Ellis模型是描述与时间无关的无屈服剪切稀释非牛顿流体的三参数模型。它常被用来代替幂律模型，且在实验测量中明显优于幂律模型。该类流体黏度计算公式为：

其中，μ_o ⁿ表示当前时刻低剪切应力时的黏度，μ_o ⁿ⁺¹表示下一时刻低剪切应力时的黏度，τ是剪切应力，τ_1/2是初始流体黏度μ_o＝μ_o/2时的剪切应力，α是经验常数。三者关系如图3所示。

当管束内存在单相Ellis模型类流体，此时，随着时间步长更新，不同时刻的非牛顿流体黏度均按照公式(27)更新。流量计算公式为：

当管束内水相与Ellis流体共存时，流量计算公式为：

此时，μ_eff的计算中，μ_o为当前时刻的流体黏度，其值随时间步长变化而按照公式(27)变化。模拟过程中，更改τ、α和μ_o，即可更改Ellis流体性质。

同样的，根据质量守恒定律构建类似AP＝B的矩阵方程，通过逐次超松弛迭代法(SOR)求解，同时需要根据时间步长Δt校正。驱替过程一直进行直到整个孔隙网络空间被侵入流体占满。然后重新计算并更新压力场和饱和度，进行下一步计算。整个过程保持注入流量Q不变。

由于流动过程中两相界面随时间推进，必须选定一个时间步长使该步长内每个两相界面均发生了适量的位移Δx，这里的“适量”是指必须在保证精度的前提下尽量减少运算次数，为此引入最小时间步长和修正时间步长。修正时间步长是指将渗流瞬态划分为长度相等的时间间隔Δt_k(k＝1,2,…)，当Δt_k足够小时，相邻时间步长的压力场变化可以认为是稳定、线性的，此时可以使用与单相流部分相同的方法(超松弛迭代法)对压力场进行求解。最小时间步长Δt_i(i＝1,2,…)是指在所有管道中的凹液面到达下一节点的时间中，选取最小的Δt_i作为这一步运算总体的时间步长，此时除到达下一节点的凹液面外的其它凹液面的位移为Δx_ij＝v_ij·Δt_min，由此可以求得该时间点的水力传导系数g_ij和两相在孔隙网络中的分布。显然最小时间步长法中时间步长Δt_i与压力降Δp和凹液面位置有关，因此每次迭代的步长取决于该次计算的具体情况，而不是都相等。最小时间步长的引入使得此模型可以使用尽量少的迭代次数得出模拟结果，保证精度的同时大大提高了计算效率。

本发明的模拟技术与传统黑油模型模拟方法中的IMPES(隐式计算压力，显式计算饱和度)模拟方法有类似之处，通过隐式有限差分方法计算孔隙流体的压力场分布，然后显示计算两相流体界面移动。孔隙网络模型可以通过GPU加速计算技术提升模型尺度，实现基于同一渗流物理模型下的从孔隙网络、室内岩心和物理模型实验、到井筒-单井油藏模型尺度的多尺度跨越。

实施例2

1、微CT扫描实验统计配位数

基于计算机高分辨断层扫描成像技术(MicroCT对样品进行了扫描及数字岩心3D重构分别用等效球法和最大球法进行了孔网构建，并对储层进行了结构特征分析。其中，吼道长度可由下式计算：

L＝D-R₁-R₂

式中，R₁，R₂分别为该吼道所连接两个孔隙的半径，单位为μm；D为两孔隙中心点的实际坐标距离，单位为μm。

配位数由软件自动统计得到。

2、核磁共振岩心分析测试统计孔喉半径分布

核磁共振信号强度与饱和岩样内部的流体氢核数量成正相关。测量方法为：测量一组样品定标，定标样品的孔隙度值分别为2％、4％、8％、15％、30％，在得到每个定标样品核磁不同回波时间的一系列T2图像，以定标样品的单位体积质子密度图像信号为横坐标，以定标样的孔隙度为纵坐标，对数据做曲线拟合，获得核磁孔隙度刻度关系式：

其中，a，b为拟合参数、为孔隙度、S为核磁图像信号大小、V为样品体积。然后测量目标岩心，将其测量图像信号值代入刻度关系式中，即可得到岩心孔隙度和整个孔隙度分布值。

同时，将核磁共振弛豫时间分布与常规岩石压汞孔径分布相结合，可以求得一个换算系数c(由压汞实验获取，其值具有地区经验性)值，将T2谱的横坐标乘以c即可获得孔径分布，如图5所示。

3.利用C++代码建立孔隙裂缝双重介质模型

本次建模储层孔喉半径呈对数正态分布(由核磁T2谱得)，通过Matlab拟合频率分布曲线，得到分布的平均值μ和标准偏差σ，得到对数正态分布随机函数：

其中，x>0，表示孔喉半径。

对应的C++生成代码：std::lognormal_distribution<double>distribution(μ,σ)。

本实施例采用SC网络模型构建孔隙网络。在SC模型网络中，每个节点连接6个相邻节点相连，因此SC网络模型的最大配位数z_max＝6。网络模型的连通概率p＝z/z_max，，使1-p倍数量的节点相互间连接的通道半径设为0，通过这种方法将部分管束去除可以生成部分连通的不规则孔隙网络模型。例如：若网络配位数z为3，则是通过随机将二维正方形网络中50％的连接通道去除后生成的模型。二维正方形网络模型的渗流临界连通概率为p_c＝66.66％，对应的平均配位数为4，本文构建了节点数为40000(节点数为200×200)的方形模型，并设置为圆形边界。通过上述对数正态分布随机分配网络中各孔道半径，网络模型的平均孔喉长度配位数6。

模拟过程中，各参数取值为：压缩系数取值为1×10^-10MPa^-1；润湿角取值为180°；两相界面张力为20mN/m²；油相黏度为100mPa/s；宾汉模型中β取值为0.3；Ellis模型中α取值为3、τ取值为100。模型为中心注入，四周采出。

4、图5为SC网络模型图，图6和图7为非稳态非牛顿流体两相驱替剖面图，图6为水驱宾汉流体流量图，图7为水驱Ellis流体流量图，中间为驱替相侵占的流动通道。

以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本发明的保护范围。

Claims (4)

1.基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其特征在于：包括以下步骤：

(1)、由核磁共振实验与恒速压汞实验获取锥状孔喉半径分布；

(2)、通过微CT扫描得到锥状孔喉长度和配位数；

(3)、根据锥状孔喉半径分布、孔喉长度和配位数，建立孔隙网络模型；

(4)、在孔隙网络模型中引入非稳态渗流理论，同时根据孔隙内流体流动、界面移动和孔隙压力扩散过程，研究孔隙型介质内部水驱油两相非稳态渗流，引用非牛顿流体，形成非牛顿两相流体驱替模拟方法，将动态网络模拟算法与非稳态渗流理论相结合模拟非牛顿两相流体驱替过程；

步骤(4)中，非稳态非牛顿两相流体模拟方法为：

若孔隙网络模型孔喉形状为圆管状，孔隙网络模型中，管束内两相流体共存时非润湿流体与润湿流体界面引起的毛细管压力p_cij通过Young-Laplace方程计算：

p_cij＝-2γcosθ/r_ij；

其中，r_ij为管束两端节点i和节点j之间管束中两相界面处的半径，γ为两相界面张力，θ为润湿角；

孔隙网络模型孔喉形状为锥状时，毛管压力为弯液面在管束中的位置的函数，称为动态毛管压力p_cij，因此，将Young-Laplace改为：

其中X_ij是与凹液面位置有关的无量纲数，即凹液面所在位置横坐标除以整个孔隙网络的长度，当管道中只有单相流体时，p_cij＝0；且管束内两相共存时，当弯液面处于管束末端时，p_cij＝0，在管束中间时具有最大值p_cij＝4γ/r_ij；

孔隙网络模型中，当存在两相弯液面时，即两相共存时，管内流量计算公式为：

其中r_ij为管束半径，p_i和p_j分别为管束两端节点i和节点j的孔隙压力，μ_eff为管束内有效黏度，其计算方法为μ_eff＝μ_IX_ij+μ_D(1-X_ij)，其中μ_I为非润湿相流体粘度，μ_D为润湿相流体黏度，初始孔隙网络模型中全部充填油相，作为被驱替相采出，模型润湿角取180，水相为非润湿相作为驱替相注入，当管束内只有非润湿相流体时，μ_eff＝μ_I，p_cij＝0，当管束内只有润湿相流体时，μ_eff＝μ_D，p_cij＝0；

非稳态非牛顿两相流体模拟中，对非牛顿流体中的宾汉模型进行研究：

(1)宾汉流体；当剪切应力超过某一静态剪切应力时，流体发生流动，其本构方程为：

τ＝τ₀+μ_bγ；

其中，τ₀为屈服应力，μ_b为宾汉流体黏度；

当管束内只有宾汉流体时，流量计算公式为：

其中r_ij0为岩心物理模型半径，在模型中，假设模型中所有管束的r_ij/r_ij0均相同，则对宾汉流体来说，上式可写为：

其中β表示与非牛顿流体流变性有关的参数，在模拟过程中，被驱替相只改变了黏度，每次迭代计算过程中，随压力变化时，黏度μ_b通过计算公式μ_bnew＝(μ_b/1-β)更新；

当管束内水与宾汉流体两相共存时，μ_eff的计算中，μ_D＝μ_b；若β＝0，则被驱替相为牛顿流体；模拟过程中，更改β，即可更改宾汉流体性质；

非稳态非牛顿两相流体模拟中，对非牛顿流体中的Ellis模型进行研究：

每次迭代计算过程中，随压力变化时，Ellis流体的黏度μ_e通过如下计算公式更新：

其中，μ_o是低剪切应力时的黏度，τ_1/2是μ_o＝μ_o/2时的剪切应力，α是经验常数；

当管束内只有宾汉流体时，流量计算公式为：

当管束内水与Ellis流体共存时，流量计算公式为：

此时μ_D＝μ_e；模拟过程中，更改τ、α和μ_o，即可更改Ellis流体性质；

根据质量守恒定律构建AP＝B的矩阵方程，通过逐次超松弛迭代法求解，同时需要根据时间步长Δt校正；驱替过程一直进行直到整个孔隙网络空间被侵入流体占满；然后重新计算并更新压力场和饱和度，进行下一步计算；整个过程保持注入流量Q不变。

2.根据权利要求1所述的基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其特征在于：步骤(1)中，建立SC模型，并在SC模型的基础上构建无序孔隙网络模型；通过核磁共振实验得到岩心孔径分布T2谱，结合与恒速压汞实验，并根据预设定量关系将所述T2谱转换得到孔喉半径分布频率，根据所述孔喉半径分布频率拟合得孔喉分布函数。

3.根据权利要求1所述的基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其特征在于：步骤(2)中，微CT扫描得到配位数与喉道长度包括以下步骤：

(2.1)利用微焦点射线源发射的锥束X射线穿透样品后投影到探测器上，同时让样品和射线源及探测器进行360°的相对旋转，采集上千帧角度的数据，然后利用计算机断层扫描成像重构方法进行3D重构，从而得到样品内外结构的高分辨3D数据及影像，从而完成利用微CT扫描实验获取CT实验数据；

(2.2)进行CT数据分析，通过图像不同灰度进行物质区分，灰度值低的区域代表物质密度低，参考灰度曲线中孔隙的灰度值进行阈值划分，在图像中将孔隙进行单独的提取分离；

(2.3)利用最大球算法，区分数字岩心三维图像中的孔隙、喉道所占空间及连通性，提取相应的孔隙、喉道结构网络模型，同时运用数理统计方法可实现对孔喉长度、孔喉体积、孔喉比、配位数、形状因子的定量提取，得到研究岩石孔喉表征的参数；

其中，喉道长度l可由下式计算：

l＝D-R₁-R₂；

式中，R₁，R₂分别为该喉道所连接两个孔隙的半径；D为两孔隙中心点的实际坐标距离。

4.根据权利要求1所述的基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法，其特征在于：步骤(3)中，建立孔隙网络模型的具体步骤如下：

(3.1)确定网络模型的大小与节点数，构建一个X×Y×Z的三维简单立方体网格；每个节点代表一个孔隙，节点与节点之间由喉道相连；由此建立起来的网络中每个代表孔隙的节点周围都有六个喉道相连，每个喉道两端与两个孔隙相连；其x方向、y方向和z方向各节点之间的间隔距离为喉道长度l，各方向节点数为nx,ny,nz，模型的边长为L_x＝(nx-1)×l、L_y＝(ny-1)×l、L_z＝(nz-1)×l；

(3.2)计算出网络模型中各节点的坐标；计算式为：

(x，y，z)＝[(i-1)l，(j-1)l，(k-1)l]；

式中i、j、k分别为x方向、y方向和z方向的节点序号；

(3.3)通过中心点位移中心节点坐标生成随机网络；

按照如下式移动各节点坐标(x，y，z)，实现各节点坐标在半径为0.5l的球面区域内随机移动，rand()为随机函数；

(3.4)在程序中设定概率为p的概率函数；

p＝z/z_max，z为配位数，表示孔喉间的连接概率，z_max为SC模型具有的最大配位数，其值为6；

如当连通概率p为66.7％时，此时对应配位数为4，在rand()函数生成的1～100整数中有66.7％的概率小于50，有另外33.3％的概率大于50；因此，将小于50的数值的管束与节点连接，并将其索引值赋为1，不满足条件的赋值为0，从而为每个管束赋予半径属性；

(3.5)为模型赋予孔喉半径属性与喉道长度属性，首先利用C++语言和科学计算库Eigen，生成2个与孔喉数量相等长度的单位向量，取名为半径向量与喉道长度向量，此时半径向量与喉道长度向量中的所有元素值均为1，表示模型中管束的半径值与喉道长度均为单位1，再通过获取的孔喉分布函数为半径向量中的每个元素赋值，这样就产生了符合孔喉半径大小的“半径向量”，完成孔喉半径赋值，从而生成一个符合岩心孔喉半径频率分布的三维张量数据体；通过对喉道长度向量中每个元素赋喉道长度值；

若孔隙网络模型孔喉形状为圆管状，则圆管状泊肃叶定律表达形式为：

孔隙网络模型孔喉形状为锥状时，则锥状孔喉的泊肃叶公式推导过程如下：

设流体通过管道微元体，通过长度为dx时的流体流动状态是层流的，并且遵循雷诺方程，此时泊肃叶定律表示为：

其中p(x)是在x，且平滑函数r(x)＝r_t+x(r_p-r_t)/l_c处的压力，l_c为锥状管道长度，r_p为孔隙半径，r_t为微元dx所在部分的喉道半径，因此：

因此可以得到锥状孔喉泊肃叶公式为：

可写为：

其中：

r_pt为孔喉比，当孔喉为均匀的管柱，孔喉比为1，α＝1，此时的锥形管束泊肃叶公式等于圆形管柱泊肃叶公式。