一种采用数字岩心模拟计算岩心渗透率的方法

摘要

本发明涉及一种采用数字岩心模拟计算岩心渗透率的方法，其包括如下步骤：a.采用波尔兹曼网格法模拟岩心多孔介质的物理模型；b.对该岩心多孔介质的物理模型进行物理场初始化及边界处理，该边界处理采用镜像边界作为波尔兹曼网格法的边界条件，以获得模拟区域及与该模拟区域具有相同构型的若干镜像区域；c.在该模拟区域内使用质量力驱使流体流动；d.根据处理后的该岩心多孔介质的物理模型，计算该模拟区域内的密度分布函数；及e. 根据该质量力和该密度分布函数模拟计算该岩心渗透率。上述方法不仅提高了岩心渗透率计算的准确度，还缩短了该岩心渗透率的计算时间。

说明书

技术领域

本发明涉及一种采用数字岩心模拟计算岩心渗透率的方法。

背景技术

岩石的孔隙度、渗透率和含油饱和度在油气储量评估中至关重要，其数据准确性会直接影响到油气资源的储量评估，甚至油气公司的股价与金融借贷能力等。

目前，岩心特征参数的传统测试方法是通过实验室直接测量获得，虽然测试的结果准确度高，但是存在测试费用高、耗时长等问题。尤其是对于低渗透油藏，准确测量岩心渗透率可能需要几天甚至数月。

数字岩心模拟是当前石油天然气勘探开发领域的最新前沿技术之一。目前采用数字岩心模拟计算岩石渗透率主要有三种方法：孔隙网格法、有限元法和波尔兹曼网格法（Lattice Boltzmann Method, LBM）。

采用孔隙网络法进行岩石渗透率模拟（Evaluation Of The Reliability Of Prediction Of Petrophysical Data Through Imagery And Pore Network Modelling, PETROPHYSICS, 50(4), 322-334,2009），其渗透率测量的准确性很大程上受制于岩石孔隙结构的简化。

有限元法对流体在岩石三维多孔介质几何模型进行数值模拟，即对流体在孔道中微观流动的不可压缩方程N-S方程进行求解，得到流体流动的速度场和压力分布场，从而可计算得到岩心渗透率（一种低渗透储层原油边界层的模拟方法，申请号：201410665585.9），虽然其模拟计算速度快，但其模拟方法复杂。

采用波尔兹曼网格法进行岩石渗透率计算，通常设置流入和流出的密度分布函数差来产生压降以驱使流体流动，并使用循环边界条件（Pressure boundary condition of the lattice Boltzmann method for fully developed periodic flows. J Zhang，DY Kwok, Physical Review E, 2006, 73(2)），然模拟体积小于岩石实际被扫描体积，计算结果对于被测岩石所在油气储层渗透率的代表性有限。综上，现有的岩石渗透率模拟计算方法存在计算准确性不高、计算耗时长、模拟方法复杂以及计算结果代表性有限等问题。

发明内容

鉴于以上内容，有必要提供一种准确、快速的采用数字岩心模拟计算岩心渗透率的方法。

一种采用数字岩心模拟计算岩心渗透率的方法，其包括如下步骤：

a.采用波尔兹曼网格法模拟岩心多孔介质的物理模型；

b.对该岩心多孔介质的物理模型进行物理场初始化及边界处理，该边界处理采用镜像边界作为波尔兹曼网格法的边界条件，以获得模拟区域及与该模拟区域具有相同构型的若干镜像区域；

c.在该模拟区域内使用质量力驱使流体流动；

d.根据处理后的该岩心多孔介质的物理模型，计算该模拟区域内的密度分布函数；及

e. 根据该质量力和该密度分布函数模拟计算该岩心渗透率。

在模拟岩心多孔介质的物理模型之前还包括扫描岩心样品，并获取岩心多孔介质的横截面的二维图像，以及将该岩心多孔介质的横截面的二维图像转换为波尔兹曼网格参数，并根据该波尔兹曼网格参数建立该岩心多孔介质的物理模型。

该岩心渗透率的计算公式为：，其中，k表示为渗透率，F表示为质量力，f_i表示为密度分布函数，e_i表示为速度函数，为弛豫时间，∆x表示为扫描分辨率。

该物理场初始化包括对该密度函数、该质量力和该弛豫时间进行初始化。

该质量力的范围为0.001-0.01N。

该模拟区域内设有i个网格，每一网格具有不同的密度分布函数，将每一网格的密度分布函数进行计算并叠加，以得到该模拟区域内的流体流速。

该模拟区域内的第i网格的密度分布函数是通过相邻网格的密度分布函数计算得到，并表示为f_i=f(f_i-1,>i+1)，其中，n≧>

在流动方向上，使用该模拟区域内的第n-1网格（或第2网格）的密度分布函数计算该镜像区域内的未知的第n+1网格（或第0网格）的密度分布函数，垂直该流体流动的方向上，使用已知的该模拟区域内的第n-1网格（或第2网格）的密度分布函数代替该镜像区域内的未知的第n+1网格（或第0网格）的密度分布函数，其中，该模拟区域内的第n-1（或第2网格）网格与该镜像区域内的第n+1（或第0网格）网格呈镜像分布。

该采用数字岩心模拟计算岩心渗透率的方法，还包括将该模拟区域内的每一网格的密度分布函数进行流速计算并叠加，以得到该模拟区域内的流速，并用于该岩心渗透率的计算。

与该模拟区域边界相邻的镜像区域和该模拟区域呈镜像分布。

相较现有技术，上述采用数字岩心模拟计算岩心渗透率的方法，通过使用镜像边界作为波尔兹曼网格法的边界条件，以获得与模拟区域具有相同构型的镜像区域，并根据该模拟区域内的密度分布函数和质量力模拟计算，以求得岩心渗透率。本发明采用镜像边界为边界条件，其能够保持流体流动始终畅通。此外，采用质量力驱动流体流动，能够使流体流动快速达到平衡。因此，本发明的采用数字岩心模拟计算岩心渗透率的方法不仅能够提高岩心渗透率计算的准确率，还能够缩短模拟计算的时间。

附图说明

图1是本发明的一较佳实施例的采用模拟岩心模拟计算岩心渗透率的流程图。

图2是现有的波尔兹曼网格法采用循环边界条件的示意图。

图3是本发明的一较佳实施例的波尔兹曼网格法采用镜像边界条件的示意图。

图4是图2中的网格采用波尔兹曼网格法设置压降驱使流动的示意图。

图5是图3中的网格采用波尔兹曼网格法设置质量力驱使流动的示意图。

主要元件符号说明

步骤S10，S20，S30，S40，S50，S60，S70

如下具体实施方式将结合上述附图进一步说明本发明。

具体实施方式

请参阅图1，本发明的一较佳实施例的采用模拟岩心模拟计算岩心渗透率的方法，其包括以下步骤：

S10、扫描岩心样品，并获得岩心多孔介质的横截面的二维图像。

在本实施例中，该扫描采用微计算机断层扫描技术（micro computed tomography，Micro-CT）进行扫描。

可以理解的，该Micro-CT扫描的分辨率的范围为0.5微米至3微米，该扫描后的岩心样品选取有效物理体积的范围为1毫米x 1毫米x 1毫米至1毫米x 1毫米x 3毫米。

在本实施例中，该Micro-CT扫描的分辨率为3微米，该扫描后的岩心样品选取有效物理体积为1毫米x 1毫米x 3毫米。

本领域技术人员能够理解的，该岩心样品的Micro-CT扫描图像是由该岩心样品图像和黑色背景组成的灰度图像，由于该岩心样品图像的灰度值与背景灰度值有明显差异。因此，对Micro-CT扫描图像进行二值化处理，可以得到相应的二值图像，并且可以进一步读取该岩心样品的横截面的二维图像。

S20、三维重构该岩心样品的横截面的二维图像，并读取表征该岩心多孔介质的参数。

可以理解的，三维重构该岩心样品的横截面的二维图像是为了获得该岩心多孔介质的三维结构的数字文件。

进一步的，图像识别并分辨该岩心样品的三维结构的数字文件，以读取表征该岩心多孔介质的参数。

表征该岩心多孔介质的参数包括基质和孔隙度。

本领域技术人员能够理解的，该基质例如是，但不局限于中孔中渗岩心。该基质还可选用低孔低渗岩心、中孔低渗岩心、高孔高渗岩心等。

在本实施例中，该基质为中孔中渗岩心，该孔隙度为11.8%。

可以理解的，三维重构该岩心多孔介质的横截面的二维图像是为了将所述该岩心多孔介质的横截面的二维图像转换为波尔兹曼网格参数，并根据该波尔兹曼网格参数建立该岩心多孔介质的物理模型。该表征岩心多孔介质的参数也即表示该波尔兹曼网格参数。

S30、将表征该岩心多孔介质的参数设置为初始条件，并根据该初始条件采用波尔兹曼网格法模拟该岩心多孔介质的物理模型。

S40、对该岩心多孔介质的物理模型进行物理场初始化及边界处理，该边界处理采用镜像边界作为波尔兹曼网格法的边界条件，以获取模拟区域及与该模拟区域具有相同构型的若干镜像区域。

该物理场初始化包括对该密度函数、该质量力和该弛豫时间进行初始化。

在本实施例中，该模拟区域是真实模拟的，该镜像区域是虚拟的。

可以理解的，相邻的镜像区域呈镜像分布，也即相邻两镜像区域呈对称分布。

如图3所示，在本实施例中，该若干镜像区域包括与该模拟区域呈镜像分布的若干虚拟的第一镜像区域及与该第一镜像区域呈镜像分布的虚拟的第二镜像区域。

可以理解的，与该模拟区域边界相邻的镜像区域与该模拟区域呈镜像分布。

该模拟区域内设有i个网格，每一网格具有不同的密度分布函数。

可以理解的，该模拟区域内的第i网格的密度分布函数是通过相邻网格的密度分布函数计算得到，可简单表示为f_i=f(f_i-1,>i+1)，其中，n≧>

可以理解的，该相邻网格包括直接与该网格边界相邻的网格及与该网格斜对角对接的网格。该网格是三维立体的（图中未示）。

将该模拟区域内的每一网格密度分布函数进行计算并叠加，以得到该模拟区域内的流体流速，并用于计算该岩心的渗透率。

可以理解的，将该模拟区域内的密度分布函数进行流速计算，以得到该模拟区域内的流速，并用于计算该岩心的渗透率。

每一镜像区域内设有与该模拟区域相对应的网格，每一网格具有不同的密度分布函数。

如图2所示，现有的波尔兹曼网格法的边界条件采用循环边界条件，其流体流动不连续。

进一步的，如图3所示，本发明采用镜像边界条件作为波尔兹曼网格法的边界条件。因此，可以始终保持流体流动的畅通性，从而提高了渗透率计算的准确度。

可以理解的，该流体密度分布函数不限于本发明图3所示的呈斜对角分布，该方向可以是水平分布、竖直分布或与水平边界存在倾斜角度分布（例如，该倾斜角度为45度，也即斜对角分布）。

S50、在该模拟区域内使用质量力驱使流体流动。

如图4所示，传统的波尔兹曼网格法采用流入和流出压降来驱使流体流动，其在流体流动过程中会形成一定的压力差分布。然而，通过流入和流出压降来驱使流体流动，其流体流动达到平衡需要较长的时间，且使用循环边界条件及压力差模拟计算岩心渗透率的准确性影响较大。

可以理解的，该流体流动达到平衡是指该模拟区域内流速不随时间变化。

可以理解的，采用该质量力而非流入和流出压降来驱使流体流动，其能够使流体流动快速达到平衡。如图5所示，本发明的各网格均设置质量力来驱使流体流动，因此流体流动能够快速达到平衡，从而缩短了计算渗透率所需的时间。

该质量力的范围为0.001-0.01N。

S60、根据处理后的该岩心多孔介质的物理模型，计算该模拟区域的密度分布函数。

该计算该模拟区域内的密度分布函数的步骤包括，将该模拟区域内的每一网格密度分布函数进行计算并叠加，以得到该模拟区域内的流体流速，并用于计算该岩心的渗透率。

进一步的，该计算该模拟区域内的流速的步骤还包括，将该模拟区域内的密度分布函数进行流速计算∑f_ie_i/∑f_i并叠加，以得到该模拟区域内的流速，并用于计算该岩心的渗透率。

可以理解的，该模拟区域内的第i网格的密度分布函数是通过相邻网格的密度分布函数计算得到，可简单表示为f_i=f(f_i-1,>i+1)，其中，n≧>

在本实施例中，在流动方向上，使用该模拟区域内的第n-1网格（或第2网格）的密度分布函数计算该镜像区域内的未知的第n+1网格（或第0网格）的密度分布函数，垂直该流体流动的方向上，使用已知的该模拟区域内的第n-1网格（或第2网格）的密度分布函数代替该镜像区域内的未知的第n+1网格（或第0网格）的密度分布函数，其中，该模拟区域内的第n-1（或第2网格）网格与该镜像区域内的第n+1（或第0网格）网格呈镜像分布。

可以理解的，该镜像区域内的第0网格与该模拟区域内的第1网格相邻，该镜像区域内的第n+1网格与该模拟区域内的第n网格相邻。

S70、根据该质量力和该密度分布函数模拟该岩心渗透率的计算公式，根据该计算公式求得该岩心渗透率。

该岩心渗透率的计算公式为：，其中，k表示为渗透率，F表示为质量力，f_i表示为密度分布函数，e_i表示为速度函数，为弛豫时间，∆x表示为扫描分辨率。

可以理解的，该速度函数e_i表示为LBM的速度函数。

本发明采用数字岩心模拟计算岩心渗透率的方法，通过使用镜像边界作为波尔兹曼网格法的边界条件，以获得若干模型对称的虚拟模拟区域，并根据该模拟区域内的密度分布函数和质量力模拟计算，以求得岩心渗透率。本发明采用镜像边界为边界条件，其能够保持流体流动始终畅通。此外，采用质量力驱动流体流动，能够使流体流动快速达到平衡。因此，本发明的采用数字岩心模拟计算岩心渗透率的方法不仅能够提高岩心渗透率计算的准确率，还能够缩短模拟计算的时间。此外，本发明可以在岩心无损的情况下，基于Micro-CT扫描可以准确分析出岩心的绝对渗透率，并动态观察流体的流动过程。因此，可以广泛运用用到煤层气、致密油等非常规油气的勘探开发。

下面通过具体的实施例及上述附图对本发明作进一步的详细描述，以下实施例可以使本专业技术人员更全面的理解本发明，但不以任何方式限制本发明。

在对比例1-2及本发明的实施例中，采用Micro CT扫描该岩心样品，其中，该Micro-CT扫描的分辨率为0.5微米，该岩心样品的有效物理体积1毫米立方，孔隙度为11.8%。

对比例1

采用传统的实验室气渗透率测量法测量，当注入气达到稳态后，测量流入/流出压力和流量，并使用达西方程计算渗透率，其测得岩心渗透率为121mD。

对比例2

使用传统的波尔兹曼网格方法，使用循环边界作为LBM的边界条件，并使用压力差为0.4kg/m³驱动流体流动，并使用岩心渗透率计算公式：，其中，k表示为渗透率，∆P表示为压力差，f_i表示为密度分布函数，e_i表示为速度函数，为弛豫时间，∆x表示为扫描分辨率。使用传统的波尔兹曼网格方法，其运算岩心渗透率的时间约为7个小时，计算岩心渗透率为89>

实施例

采用本发明采用数字岩心模拟计算岩心渗透率的方法，使用镜像边界边界作为LBM的边界条件，并使用质量力为0.01N驱动流体流动，并使用该岩心渗透率的计算公式：，其中，k表示为渗透率，F表示为质量力，f_i表示为密度分布函数，e_i表示为速度函数，为弛豫时间，∆x表示为扫描分辨率。本发明采用数字岩心模拟计算岩心渗透率的方法，其运算岩心渗透率的时间时间约为3个小时，计算岩心渗透率为132>

性能测试

将对比文件1-2及实施例模拟计算该岩心样品的渗透率、运算该岩心样品的渗透率的时间进行测试，并比采用传统LBM模拟计算的渗透率及与的与实验室实测的渗透率进行比较，具体的测试结果如表一所示。

表一

如表一所示，对比例2计算的岩心渗透率为89mD，远小于对比例1测得的岩心渗透率，其相对误差为26.4%。然本发明实施例1计算的岩心渗透率为132mD，略高于对比例1测得的岩心渗透率，其相对误差为9.1%。因此，测试结果说明利用本发明的采用数字岩心模拟计算岩心渗透率的方法计算岩心渗透率，其计算结果更接近实际测量的岩心渗透率。相较于对比例2，本发明的特征岩心渗透率模拟计算结果的误差可减小至10%以内，且计算速度提高1倍以上。

本发明采用数字岩心模拟计算岩心渗透率的方法，通过使用镜像边界作为波尔兹曼网格法的边界条件，以获得若干与模拟区域具有相同构型的虚拟的模拟区域，并根据真实的模拟区域内的密度分布函数和质量力模拟计算，以求得岩心渗透率。不仅能够提高岩心渗透率测量的准确度，还缩短了模拟计算的时间。此外，与现有的数字岩心模拟计算方法相比，在不改变被扫描岩心物理体积的前提下，本发明采用的数字岩心模拟可扩大表征体元（ Representative Elementary Volume， REV）体积。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，以上实施方式仅是用于解释权利要求书。然本发明的保护范围并不局限于说明书。任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或者替换，都包含在本发明的保护范围之内。