CN117057271B

CN117057271B - 一种基于vof的多相流体渗流过程模拟方法

Info

Publication number: CN117057271B
Application number: CN202311027265.6A
Authority: CN
Inventors: 唐雁冰; 李星甫; 周利华; 杨培强; 李闽; 吴飞; 卢丙; 杨鑫; 张柏宇; 杜翔宇
Original assignee: Southwest Petroleum University
Current assignee: Southwest Petroleum University
Priority date: 2023-08-15
Filing date: 2023-08-15
Publication date: 2024-03-01
Anticipated expiration: 2043-08-15
Also published as: CN117057271A

Abstract

本申请提供了一种基于VOF的多相流体渗流过程模拟方法，涉及油气开发技术领域，该方法包括：S1:获取岩心样品的实际孔隙度φ和实际渗透率k；S2:利用扫描技术获取岩心样本在不同断面上的二维图像数据；S3:根据二维图像数据，获得岩心三维数据体；S4:根据实际孔隙度φ与二维图像数据的关系，获得岩心样品的模拟孔隙度；根据实际渗透率k与所述实际孔隙度φ、二维图像数据的关系，获取岩心样品的模拟渗透率；S5:利用所述模拟孔隙度和所述模拟孔隙度得到岩心渗透率三维数据体；S6:根据岩心渗透率三维数据体赋值三维立方体模型，得到高精度数字岩心模型。本申请提供的方法适用于低渗、超低渗储层难以获得孔喉数据的技术问题。

Description

一种基于VOF的多相流体渗流过程模拟方法

技术领域

本申请涉及油气开发技术领域，尤其是涉及一种基于VOF的多相流体渗流过程模拟方法。

背景技术

在油气田开发及渗流力学领域，常规油气藏数值模拟技术主要用于研究油气开采过程中地下油气的分布和运动状况，但就实际应用而言，目前存在以下问题：

在多相渗流中，相界面的存在给数值模拟技术的应用带来很大障碍。首先，界面将渗流区域分开，单相流体渗流的数值计算技术难于直接应用。再者，界面的形状有时较为复杂，很多不稳定性机制控制着界面的运动特性，由于相界面及其随时间的迁移、变形、破碎与融合等复杂变化，使两相流的渗流特性比单相流复杂得多，相界面的迁移特性是两相流研究中的一个重要问题(高一娟.微通道内两相流界面追踪的数值模拟[D].天津大学,2009.)。对于多相流体的研究，常规油气藏数值模拟技术计算通常比较难，尤其是无法准确进行界面追踪。

发明内容

本申请的目的在于提供一种基于VOF的多相流体渗流过程模拟方法，以解决如何更准确获取岩心尺度的两相多组分渗流过程，帮助分析两相驱替过程的分布状态、动用状况以及渗流过程影响因素的技术问题。

本申请的另一目的在于提供一种基于VOF的多相渗流模型构建系统。

为了上述目的，本申请提供以下技术方案：

第一方面，本申请提供一种基于VOF的多相流体渗流过程模拟方法，包括：

S1:利用岩心样品和扫描技术，获取所述岩心样品的孔隙度三维数据体和渗透率三维数据体；

S2:建立所述岩心样品的三维立方体模型，根据所述孔隙度三维数据体和渗透率三维数据体赋值所述三维立方体模型，得到所述岩心样品的高精度数字岩心模型；

S3:根据动量守恒方程、质量守恒方程以及相方程，构建多相渗流数学模型；

S4:根据所述岩心样品的高精度数字岩心模型求解多相渗流数学模型，得到所述岩心样品的室内数字化多相渗流模拟结果。

进一步地，在本申请的一些实施例中，所述动量守恒方程为：

式中，G_m为传导率；ρ_o、ρ_g以及ρ_w分别为油、气、水三相的密度，单位为kg/m³；μ_o、μ_g、μ_w分别为油、气、水的粘度，单位Pa·s；k_o、k_g、k_w分别为油、气、水三相流体的渗透率，单位为m²；p_co、p_cg、p_cw分别为油、气、水的毛管压力；和/>为油、气、水三相的真实渗流速度。

进一步地，在本申请的一些实施例中，所述质量守恒方程为：

式中，t为时间，单位为s；ρ_o、ρ_g、ρ_w为油、气、水三相流体的密度，单位为kg/m³；V_o、V_g、V_w分别表示微元体里面油、气、水三相的体积，单位为m³；和/>为油、气、水三相的真实渗流速度，单位为m/s。

进一步地，在本申请的一些实施例中，所述相方程为：

式中，S_o、S_g和S_w分别为油、气、水三相的饱和度，和/>为油、气、水三相的真实渗流速度，单位为m/s；t为时间，单位为s。

进一步地，在本申请的一些实施例中，所述求解包括基于多相渗流数学模型获取所述岩心样品在渗流实验中各相流体的压力场和饱和度场。

进一步地，在本申请的一些实施例中，所述孔隙度三维数据体和渗透率三维数据体的获得包括以下步骤：

S11:获取待模拟储层的岩心样品，并利用所述岩心样本获取所述岩心样品的实际孔隙度φ和实际渗透率k；

S12:利用扫描技术获取所述岩心样本在不同断面上的二维图像，获得二维图像数据；

S13:根据多个所述断面上的二维图像数据，获得岩心三维数据体；

S14:根据岩心的实际孔隙度φ与二维图像数据的关系，获得所述岩心样品的孔隙度转化系数；根据所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系，获取所述岩心样品的渗透率转化系数；

S15:在所述岩心三维数据体上应用所述孔隙度转化系数和所述渗透率转化系数，以及实际孔隙度与二维图像数据的关系、实际渗透率与实际孔隙度及二维图像数据的关系，得到岩心孔隙度三维数据体和岩心渗透率三维数据体。

进一步地，在本申请的一些实施例中，所述扫描技术为核磁共振扫描技术或者CT扫描技术中的任一种。

进一步地，在本申请的一些实施例中，当所述扫描技术为核磁共振扫描技术时，所述二维图像数据包括图像数据和断面坐标；所述实际孔隙度φ与二维图像数据的关系为φ＝a_MRI×V_MRI；其中，V_MRI为图像数据，α_MRI为转换系数；所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系：k＝β_MRI×V_MRI ²×φ；β_MRI为转换系数；

当所述扫描技术为CT扫描技术时，所述二维图像数据包括图像数据和断面坐标；所述实际孔隙度φ与二维图像数据的关系为φ＝α_CT×1/V_CT；其中，V_CT为图像数据，α_CT为转换系数；所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系：k＝β_CT×(1/V_CT ²)×φ；β_CT为转换系数。

进一步地，在本申请的一些实施例中，所述图像数据为岩石密度或者灰度或者像素数中的任一种。

第二方面，本申请提供一种基于VOF的多相流体渗流过程模拟方法在油气藏勘探和开发技术领域中的应用。

本申请提供一种基于VOF的多相流体渗流过程模拟方法，在多相渗流数学模型中引入相体积分数，并以相体积分数作为变量来对计算域内相间界面的追踪、描述流动相的传输行为，可以更加准确的模拟得到两相流体界面，从而能更精确获取岩心尺度的两相多组分渗流过程，有助于分析两相驱替过程的分布状态、动用状况以及渗流过程的影响因素。此外，本申请提供的模拟方法基于高精度的数字岩心模型，该数字岩心模型利用扫描技术如核磁共振成像MRI扫描技术或者CT扫描技术获取岩石样品的断面二维图像数据如图像像素、灰度，用以表示岩石样品的岩石密度；同时直接根据岩石密度与孔隙度、渗透率的关系，直接获取岩石样品的被扫描处的孔隙度、渗透率等参数，解决了现有技术中无法直接获取岩心样品油、气、水三相相对渗透率等物理性质的缺陷，且该方法无需借助孔喉长度来获得孔隙度、渗透率，其构建过程简单，尤其适用于低渗储层(渗透率10mD～100mD、孔隙度15％～20％)、超低渗储层(渗透率低于10mD、孔隙度不高于15％)获取的岩心样品的孔隙度和渗透率的直接获得，解决了现有技术中低渗储层、超低渗储层的岩心样品的孔喉性能获取难的缺陷。本申请将将岩心的数字化孔隙网络模型这一数字岩心模型与多相渗流数值模拟方法相结合，可模拟再现室内岩心流动(驱替、自吸等)真实实验过程和结果，且数值模拟能跟岩心流动实验的图像分析结果进行直接对比和分析，确保模拟结果的真实可靠。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。

图1为本申请提供的基于VOF的多相流体渗流过程模拟方法的流程示意图；

图2为本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例1得到的岩心样品JY-27水驱前后整体核磁成像图(依次为饱和、0.5PV、1.0PV、1.5PV、2.0PV)；

图3本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例1得到的岩心JY-27饱和、驱替后核磁扫描T₂谱；

图4本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例1得到的岩心样品JY-27在不同驱替速度下油水分布特征示意图；

图5本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例2得到的岩心样品ZJ-10水驱前后整体核磁成像图(左为驱替前，右为驱替后)；

图6本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例2得到的岩心样品ZJ-10核磁分层成像图(上排为驱替前、下排为水驱后，左侧为入口端，右侧为出口端)；

图7本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例2得到的岩心样品ZJ-10在不同驱替速度下油水分布特征示意图；

图8本申请提供的基于VOF的多相流体渗流过程模拟方法中实施例2得到的岩心样品ZJ-10水驱油采出程度随注入PV变化关系图。

具体实施方式

下面将结合实施例对本申请方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

发明人在本申请中提出了一种基于VOF的多相流体渗流过程模拟方法，参阅图1，包括：

其中，所述动量守恒方程为：

所述质量守恒方程为：

所述相方程为：

式中，G_m为传导率；ρ_o、ρ_g以及ρ_w分别为油、气、水三相的密度，单位为kg/m³；μ_o、μ_g、μ_w分别为油、气、水的粘度，单位Pa·s；k_o、k_g、k_w分别为油、气、水三相流体的渗透率，单位为m²；p_co、p_cg、p_cw分别为油、气、水的毛管压力；和/>为油、气、水三相的真实渗流速度；式中，t为时间，单位为s；ρ_o、ρ_g、ρ_w为油、气、水三相流体的密度，单位为kg/m³；V_o、V_g、V_w分别表示微元体里面油、气、水三相的体积，单位为m³；/>和/>为油、气、水三相的真实渗流速度，单位为m/s；式中，S_o、S_g和S_w分别为油、气、水三相的饱和度。

由于岩石在渗流过程中，根据牛顿第二定律，有如下表达式：

其中，F_pres、F_grav、F_visc、F_cap分别为流体微元受到周围流体的净力、重力、粘滞力以及毛管力，其中压力作为动力，粘滞力、惯性力为阻力，而重力和毛管压力既有可能作为动力，也有可能作为阻力，所有力的单位均为N(忽略重力)。对一个微元体里面的单相渗流流体(比如油、气、水三相流体中的任意一相)进行受力分析，x方向如式(2a)所示，y方向如式(2b)所示，z方向如式(2c)所示，此时F_cap＝0，m为微元体流体的质量，其单位为kg；a_x为x方向的加速度，a_y为y方向的加速度，a_z为z方向的加速度，加速度单位为m/s²，式(2a)、式(2b)、式(2c)如下所示：

F_presx+F_viscx＝ma_x (2a)

F_presy+F_viscy＝ma_y (2b)

F_presz+F_viscz＝ma_z (2c)

其中流体微元受到周围流体在x、y、z方向上的净受力F_presx、F_presy、F_presz可以表示为如下式:

其中p为压力，单位为Pa；A为微元体的面积，单位为m²；dV为微元体的体积(dV＝dxdydz)，单位为m³；粘滞力作为阻力，单位为N；粘滞力在x、y、z方向上(F_viscx、F_viscy、F_viscz)可以表示为如下公式:

其中，μ为单相流体的粘度，单位为Pa·s；k_x、k_y、k_z分别为x、y、z方向的绝对渗透率，单位为m²；w_x、w_y、w_z分别为x、y、z方向的真实渗流速度，单位为m/s；

基于上述式子，式(2a)、(2b)、(2c)可以进一步转换为：

将式(5a)、(5b)、(5c)写为如下的通式：

其中压力梯度可以用如下公式来表示：

则式(6)可以进一步推导为：

其中真实渗流速度的物质导数/>可以表示为：

在式(8)中引入连续性方程，其表达式如下：

再在上述式(10)上述方程两边同时乘以速度则可以得到：

联立公式(8)、公式(9)、公式(11)可以得到：

基于单相流体的动量守恒方程：

忽略掉对流这一项得到

因此，油、气、水三相对应的单相流体动量守恒方程如下：

式中，φ为岩石的孔隙度，ρ_o、ρ_g和ρ_w分别为油、气和水的密度，单位为kg/m³；μ_o、μ_g、μ_w分别为油、气、水的粘度，单位Pa·s；和/>为油、气、水三相的真实渗流速度，单位为m/s；k_o、k_g、k_w分别为油、气、水三相流体通过真实渗流速度计算得到的渗透率，单位为m²。

当岩心样品在进行渗流试验时，将岩心样品分为若干发生渗流的微元体，根据质量守恒定律，通过微元体控制面净流入该物理量的通量与微元体单位时间内流体质量的增加量相等，方向是流出为正，流入为负，其中单位时间内油、气、水质量增加量的公式如下：

式中，t为时间，单位为s；ρ_o、ρ_g、ρ_w为油、气、水三相流体的密度，单位为kg/m³；其中单位时间内微元体控制面净流入通量的表达式如下：

A_o、A_g、A_w分别表示微元体里面油、气、水三相当前控制面的面积，单位为m²；其中对于面积有如下的转换公式：

dA_o＝φS_odA (18a)

dA_g＝φS_gdA (18b)

dA_w＝φS_wdA (18c)

由此质量守恒方程可以转换为：

其中，V为微元体的体积，单位为m³；A为控制面的面积，单位为m²。

因此可以得到：

对上述式(19)引入源汇项Q(表示某一节点内可能有外部流体注入，或有流体从微元体中流出；否则Q＝0，源汇项为质量流量，单位为kg/s)，可以得到式(20)：

其中ρ_m为混合流体的密度，其单位为kg/m³，其表达式如下：

ρ_m＝S_oρ_o+S_gρ_g+S_wρ_w (22)

而w_m为混合流体的渗流速度，其单位为m/s，其表达式如下：

而对于油、气、水三相流体的多相渗流情况，动量守恒方程还需要考虑毛管压力的作用，如果油相占优，那么此时三相毛管压力可以表示为：

p_c＝p_cog+p_cow (24)

如果气相占优，那么此时三相毛管压力可以表示为：

p_c＝p_cog+p_cwg (25)

如果水相占优，那么此时三相毛管压力可以表示为：

p_c＝p_cow+p_cwg 26)

p_cog、p_cow、p_cwg分别为油气、油水以及气水直接的毛管压力，其表达式如下：

其中，s_og、s_ow、s_wg分别为油气、油水、气水之间的界面张力，单位为N/m；q_og、q_ow、q_wg分别为油气、油水、气水之间的润湿接触角，无量纲；r为喉道半径，单位为m。

对于油相的动量守恒方程，在x，y，z方向上分别有:

其中w_ox、w_oy、w_oz分别为油相在x，y，z方向的真实渗流速度，单位为m/s。

对于气相的动量守恒方程，在x，y，z方向上分别有：

其中w_gx、w_gy、w_gz分别为气相在x，y，z方向的真实渗流速度，单位为m/s。

对于水相的动量守恒方程，在x，y，z方向上分别有：

其中w_gx、w_gy、w_gz分别为水相在x，y，z方向的真实渗流速度，单位为m/s。

方程两边同时乘以此时假设毛管力为阻力，以x方向为例，此时油、气、水三相公式可以表示为如下：

将上述式子相加可以得到x方向上的动量守恒方程：

其中x方向传导率G_mx可以表示为：

同时：

在y，z方向上同理可进行上述推导，即可得到针对多相渗流的动量守恒方程:

其中，任意网格内的油、气、水三相相对渗透率k_ro、k_rg和k_rw可采用修正Brooks-Corey模型进行计算，其公式如下，相对渗透率无量纲：

式中，n_o、n_g、n_w分别为油、气、水的相渗指数(取值范围为1～5)，相渗指数无单位；S_or、S_gc、S_wc分别为控制体内残余油、残余气和束缚水的饱和度值，饱和度无单位。

其中，相方程为用于求解多相渗流中各组分的体积分数的连续方程，其具体表达如下：

联立上述动量守恒方程、质量守恒方程、相方程即可得到多相渗流数学模型。

然后基于岩心样品的数字岩心模型对多相渗流数学模型进行求解，即可得到多相渗流的岩心样品室内渗流模拟结果，其求解包括压力场求解和饱和度场求解，其具体求解过程包括以下步骤：

(1)压力场求解：

对于式(22)的第一项，

/>

其中，C_o、C_g、C_w分别为油、气、水三相的压缩系数，压缩系数无量纲；这一项可以进一步简化为:

其中，t为当前时刻，t+Δt为下一时刻，单位均为s。而对于公式(35)，可以忽略掉惯性力项(即不考虑时间项)，得到如下公式：

简化后，对于公式(21)的第二项可以写为，

其中，z_c为配位数，配位数无量纲；A_ij为相邻网格的面积，单位为m²。联立公式(23)、公式(40)、公式(42)可以得到：

所有网格均可将当前时刻t和下一时刻t+Δt的不同参数代入公式(42)中，从而可得到由N个与上述公式相似的方程组成的方程组。根据方程组的下标，可将方程组整理为矩阵形式，即A^t+ΔtX^t+Δt＝B^t：

式中，A^t+Δt为N×N大小的与流体水力传导率有关的稀疏矩阵(N为模型的网格数量)，X^t+Δt和B^t是长度为N的两个向量，X^t+Δt为下一时刻t+Δt的压力场向量，B^t为与当前时刻t压力场和边界条件有关的向量，求解以上矩阵方程即得到下一时刻模型中流体的压力场。

(2)求解饱和度场：基于式(37)，引入源汇项Q_o、Q_g以及Q_w，其中Q_o、Q_g以及Q_w为油、气、水三相的体积流量，单位为m³/s，表达式如下：

进一步差分后可以得到油、气、水的分相渗流方程，用以计算网格中每一相流体的饱和度：

其中，为当前时刻，下一时刻油相的饱和度；/>为当前时刻，下一时刻气相的饱和度；/>为当前时刻，下一时刻水相的饱和度；饱和度均无量纲，此时可以得到每一相模拟求解饱和度的公式：/>

其中，需要说明的是，岩心样品的高精度数字岩心模型的构建包括以下步骤：

S15:在所述岩心三维数据体上应用所述孔隙度转化系数和所述渗透率转化系数，以及实际孔隙度与二维图像数据的关系、实际渗透率与实际孔隙度及二维图像数据的关系，得到岩心孔隙度三维数据体和岩心渗透率三维数据体；

S16:建立所述岩心样品的三维立方体模型，根据所述岩心孔隙度三维数据体和岩心渗透率三维数据体赋值所述三维立方体模型，得到所述岩心样品的高精度数字岩心模型。

需要说明的是，所述岩心样品经过洗油洗盐并且在一定温度下，如80℃下烘干后的岩心样品。其中，所述扫描技术可以为核磁共振扫描技术或者CT扫描技术中的任一种，所述岩心样品在饱和地层水之后进行断面扫描。其岩心样品进行扫描时，其断面位置和断面数可以根据扫描时所采用的仪器的精度进行调整，精度高，则其断面数可以适当提高；精度低，则其断面数可以适当降低。其中获得的二维图像数据包括图像数据，即图像的像素体数据，和断面坐标，即切片位置坐标。获得的二维图像数据，可以保存为TXT文本，以便于进行后续数据处理。

当所述扫描技术为核磁共振扫描技术时，所述实际孔隙度φ与二维图像数据的关系为φ＝α_MRI×V_MRI；其中，V_MRI为图像数据，α_MRI为转换系数；所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系：k＝β_MRI×V_MRI ²×φ；β_MRI为转换系数；

当所述扫描技术为CT扫描技术时，所述实际孔隙度φ与二维图像数据的关系为φ＝a_CT×1/V_CT；其中，V_CT为图像数据，α_CT为转换系数；所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系：k＝β_CT×(1/V_CT ²)×φ；β_CT为转换系数。

其中，所述图像数据，即像素体数据可以为岩石密度或者灰度或者像素数中的任一种。由于图像的灰度或者像素数可以用于表示岩石密度的大小，因此在本申请中直接用灰度或者像素数直接表示相对的岩石密度。

由于扫描得到的数据体难以正好落在便于处理的范围内，因此在所述步骤S3中，在根据多个所述断面上的二维图像数据，获得岩心三维数据体之前，先对所述二维图像数据进行粗化处理和/或插值处理。如，当数据体中的节点数达到上千万甚至上亿时，对所述二维数据图像进行粗化。其中，所述二维图像数据经所述粗化处理后，其网格数量不低于100万且不高于300万；而当所述二维图像数据的节点数低于100万个时，对所述二位图像数据进行插值处理，利用已知邻近像素点的灰度值(或rgb图像中的三色值)来产生未知像素点的灰度值，以便由原始图像再生出具有更高分辨率的图像。

其中，所述实际孔隙度φ和实际渗透率k为任一所述利用扫描技术扫描的岩心样品的切片的实际孔隙度φ和实际渗透率k；或者；所述实际孔隙度φ和实际渗透率k为整个岩心样品的测试得到的实际孔隙度和实际渗透率。

在一些实施例中，所述插值算法包括三线性插值算法或克里金插值算法。

在一些实施例中，其特征在于，在步骤S11中，当所述扫描技术为CT扫描技术时，所述岩心样品通过洗油、洗盐以及烘干处理后并测量实际孔隙度φ和实际渗透率k以及岩石样本的长度、直径；

当所述扫描技术为核磁共振扫描技术时，所述岩心样品通过洗油、洗盐以及烘干处理后，再进行饱和地层水处理，并测量实际孔隙度φ和实际渗透率k以及岩石样本的长度、直径。

下面将结合实施例对本申请的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

实施例1

以长庆油田低渗岩心JY-27为研究对象，将该岩心样品洗油洗盐并且在80℃烘干后，测量岩心样品的长度、直径、实际孔隙度和实际渗透率；该岩心样品的长度为5cm、直径为2.5cm、实际孔隙度φ＝11.70％、实际渗透率k＝0.1056mD。

核磁共振测量采用核磁共振分析仪(由苏州纽迈分析仪器股份有限公司提供，其型号为MacroMR12-150V-I)，参照SY/T6490-2014《岩样核磁共振参数实验室测量规范》标准，在扫描岩心过程中，将长为5cm直径为2.5cm的长庆油田标准岩心ZJ-10置于恒定的磁场内，对岩心不同部位进行切片处理，可以得到不同位置的岩心MRI二维图像。根据仪器扫描精度，在仪器可处理范围内选择合适位置和断面数(例如断面数为6)，获取岩心不同位置的核磁共振MRI图像，将切片位置坐标与像素体数据保存与TXT文本中。

由于核磁成像扫描得到的数据精度不足，需要将获取的二维图像TXT文本数据进行插值处理，插值后得到扫描区的三维数据体{V_MRIi}。

并利用与实施例1相同的步骤，利用真实孔隙度φ与核磁共振成像(<V_MRI>)的正比关系；实际渗透率k与岩心实测孔隙度φ、核磁共振成像(<V_MRI>)有如下关系式：k＝β_MRI(<V_MRI>)²φ得到岩心孔隙度三维数据体和岩心渗透率三维数据体。再将得到的岩心孔隙度三维数据体和岩心渗透率三维数据体赋值岩心样品的三维立方体模型，即可得到高精度数字岩心模型。

所述JY-27岩心样品的核磁整体成像图如图2所示.两块岩心在驱替前，模拟油密集地分布在该岩心内部，且该分布较为均匀；随着注入体积的增加，相比于岩心出口端和岩心中间部分，靠近入口端岩心的含油量降低最快，注入体积达到0.5PV之后，含油量下降已经大幅度趋缓，岩心内部的残余油向出口端集中；当注入体积达到1.5PV之后，核磁分层图像的信号已经不再随驱替时间的增加而改变，表明此时已经趋近于残余油状态。驱替结束后，核磁共振整体成像以及分层成像均显示此时该岩心内部的少许剩余油位于岩心出口端，中间部分和岩心入口端基本上不再残存剩余油。

从核磁扫描T₂谱图(图3)也可以看到，岩心JY-27中大孔小孔里面的油均被驱替了出来，且驱替效率随着注入PV数的增加(0.5PV、1.0PV、1.5PV和2.0PV)依次为22.77％、46.85％、54.54％和55.13％,表明在驱替开始的很短一段时间内，岩心内部饱和油随即被大量驱出，之后出油速率逐渐趋缓，直至岩心达到残余油状态。

模拟过程中岩心数字化模型中饱和油，将水以不同注入压力从岩心左端面中心入口处注入模型中驱替模型中的油。设置模型右端出口端的压力为大气压0.1MPa，将上述条件带入多相流体渗流数学模型中，采用本发明提出的多相渗流数学模型和数值模拟方法，进行岩心的水驱油模拟分析，岩心水驱实验与模拟对比见表1。

表1岩心JY-27水驱实验与模拟对比表

岩心样品JY-27的模拟结果如图4所示。模拟和实验结果均表明，岩心JY-27的水驱油过程中，注入前缘推不均匀，这是因为注入流体先经过有斜层理、裂缝充填的层状分布的高渗区，而该区域并未连通注入端与出口端，因此无法有效的形成优势通道。同样地，随着注入PV数增加，高渗层状裂缝区域被注入水完全充填后，流体进入均匀的低渗区，高渗区中的大量流体分布变相的增加了注入流体与岩石的接触面积，使得水更容易进入岩心JY-27置换被驱替相，采出程度较高。在驱替中后期，由于该岩心孔喉分布均匀(除裂缝高渗区外)，整体孔喉动用程度较高，而高渗区大量注入液形成的高压力降会源源不断的驱替油，从而增加波及效率。

实施例2

本实施例以岩心样品ZJ-10为示例，该岩心样品洗油洗盐并且在80℃烘干后，测量岩心样品的长度、直径、实际孔隙度和实际渗透率；该岩心样品的长度为5cm、直径为2.5cm、实际孔隙度φ＝14.35％、实际渗透率k＝1.1311mD。其依然采用核磁共振扫描技术扫描得到核磁整体成像图和核磁分层成像图，如图5、图6所示。

从图中可以看出，岩心样品ZJ-10在初始状态下模拟油密集地分布在该岩心内部，且该分布较为均匀，这点可以从核磁共振整体成像和核磁共振分层成像中均可以得出结论；水驱后岩心靠近入口端的含油量减少明显，并且越靠近入口端，含油量减小越明显；驱替结束后，核磁共振整体成像以及分层成像均显示此时岩心内部的剩余油大都位于出口端，入口端基本上不再残存剩余油，而中间部分的剩余油含量介于入口端和出口端之间。

岩心样品ZJ-10气测渗透率为1.131mD，孔隙度为14.35％，岩心样品较均质，高渗区呈点状分布。模拟过程中岩心数字化模型中饱和油，将水以不同注入压力从岩心左端面中心入口处注入模型中驱替模型中的油。设置模型右端出口端的压力为大气压0.1MPa，将上述条件带入多相流体渗流数学模型中，采用本发明提出的多相渗流数学模型和数值模拟方法，进行岩心的水驱油模拟分析，岩心样品ZJ-10水驱实验与模拟对比见表2。

表2岩心ZJ-10水驱实验与模拟对比表

岩心样品ZJ-10模拟表明，在驱替前，模拟油密集地分布在该岩心中间部分和岩心的出口段，且该分布较为均匀；水驱后岩心靠近入口端的含油量减少明显，并且越靠近入口端，含油量减小越明显；驱替结束后，两相流体分布显示此时岩心内部的剩余油大都位于出口端，入口端基本上不再残存剩余油，而中间部分的剩余油含量介于入口端和出口端之间。如图7所示，模拟该岩心的水驱油过程中，注入前缘推进较均匀，但由于该岩心高渗区程点状分布这一特征，初始注入主要均匀的从大孔通过并突破岩心出口端，导致除大孔喉外的岩心整体孔喉动用程度并不高，采出程度只有图8示的26.22％。随注入PV数增加，部分注入水虽然可以有效的进入小孔喉置换油，但这种置换效率并不高，因为大孔喉形成的优势通道会使得压力降快速降低，从而降低波及效率的增长。

综上，本申请能够较为完整地保留岩心微观物理特征，精度较高，克服了低渗、超低渗岩石孔喉数据难以获取的难题，也能够计算油、气、水三相相对渗透率等传统物理实验无法直接测量的物理性质。

最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。

Claims

1.一种基于VOF的多相流体渗流过程模拟方法，其特征在于，包括：

S4:根据所述岩心样品的高精度数字岩心模型求解多相渗流数学模型，得到所述岩心样品的室内数字化多相渗流模拟结果；

其中，所述动量守恒方程为：

式中，G_m为传导率；ρ_o、ρ_g以及ρ_w分别为油、气、水三相的密度，单位为kg/m³；μ_o、μ_g、μ_w分别为油、气、水的粘度，单位Pa·s；k_o、k_g、k_w分别为油、气、水三相流体的渗透率，单位为m²；p_co、p_cg、p_cw分别为油、气、水的毛管压力；和/>为油、气、水三相的真实渗流速度；

所述质量守恒方程为：

式中，t为时间，单位为s；ρ_o、ρ_g、ρ_w为油、气、水三相流体的密度，单位为kg/m³；V_o、V_g、V_w分别表示微元体里面油、气、水三相的体积，单位为m³；和/>为油、气、水三相的真实渗流速度，单位为m/s；

所述相方程为：

式中，S_o、S_g和S_w分别为油、气、水三相的饱和度，和/>为油、气、水三相的真实渗流速度，单位为m/s；t为时间，单位为s；

2.根据权利要求1所述的一种基于VOF的多相流体渗流过程模拟方法，其特征在于，所述求解包括基于多相渗流数学模型获取所述岩心样品在渗流实验中各相流体的压力场和饱和度场。

3.根据权利要求1所述的一种基于VOF的多相流体渗流过程模拟方法，其特征在于，所述孔隙度三维数据体和渗透率三维数据体的获得包括以下步骤：

4.根据权利要求3所述的一种基于VOF的多相流体渗流过程模拟方法，其特征在于，所述扫描技术为核磁共振扫描技术或者CT扫描技术中的任一种。

5.根据权利要求4所述的一种基于VOF的多相流体渗流过程模拟方法，其特征在于，当所述扫描技术为核磁共振扫描技术时，所述二维图像数据包括图像数据和断面坐标；所述实际孔隙度φ与二维图像数据的关系为φ＝α_MRI×V_MRI；其中，V_MRI为图像数据，α_MRI为转换系数；所述实际渗透率k与所述实际孔隙度φ、二维图像数据的关系：k＝β_MRI×V_MRI ²×φ；β_MRI为转换系数；

6.根据权利要求3所述的一种基于VOF的多相流体渗流过程模拟方法，其特征在于，所述图像数据为岩石密度或者灰度或者像素数中的任一种。