CN106442271A - 岩心渗透率模拟方法及装置 - Google Patents

Abstract

本发明提供了一种岩心渗透率模拟方法及装置，涉及油田开发的技术领域，包括以下步骤：读取所述岩心的CT图像数据；对所述CT图像数据进行二值化处理，生成所述岩心多孔介质的二值化CT图像数据；对所述二值化CT图像数据进行三维重建，生成所述岩心的三维体数据；通过所述三维体数据获取所述岩心的表征体元数据；根据所述表征体元数据建立所述岩心的物理模型，并由所述物理模型计算宏观物理量，其中，所述物理模型分布函数的存储采用稀疏矩阵存储算法；根据所述宏观物理量计算所述岩心渗透率。本发明解决了现有岩心渗透率计算效率低的问题，在提高运算速度的同时提高了渗透率的计算精度。

Description

岩心渗透率模拟方法及装置

技术领域

本发明涉及计算流体力学和油田开采技术交叉领域，尤其是涉及一种岩心渗透率模拟方法及装置。

背景技术

多孔介质的渗透率是流体在多孔介质中流动能力的重要体现，也是多孔介质流动中的基本参数。渗透率分为绝对渗透率、有效渗透率和相对渗透率，根据达西定律，绝对渗透率公式为

其中，Q是指在压差ΔP的作用下，通过岩心的流量，cm³/s；A是指垂直于流动方向的岩心截面积，cm2；L为岩心的长度，cm；μ是岩心的流体粘度，mPa·s；ΔP岩心两端的压力差，atm。K的单位通常用达西(D)表示，其中1D＝1μm²。而且，绝对渗透率是多孔介质的固有特性，与流体的性质无关。

有效渗透率是指多相流体共存流动时，岩石让其中某一相流体通过的能力。定义如下

其中Q_i表示流体i通过岩心的流量，u_i是流体i的粘性。在实际应用中，常使用无量纲化的相对渗透率，某一相流体的相对渗透率是指该相流体的有效渗透率与绝对渗透率的比值，公式为

目前研究多孔介质相对渗透率的方法主要由实验法、理论和数值模拟三大类。在数值模拟方面，多孔介质孔隙尺度的数值研究主要包括孔隙网络模型(Pore NetworkModel，简称PNM)和直接数值模拟(Direct Numerical Simulation，简称DNS)，孔隙网络模型是将复杂的多孔结构抽象为简单的几何形状，并且在计算的过程中采用了一系列的假设，譬如在准静态模型中，不考虑网络两端的粘性压降，流体的驱替形式完全由毛管力控制等等。上述的简化以及假设使得孔隙网络模型与真实的流动状态偏差较大。DNS结合CT(Computed Tomography，简称CT)扫描的数字岩心成功地克服了以上的不足，现有技术通常使用DNS结合二维CT数据来模拟岩心的绝对渗透率，由于二维截面数据不能反映岩心的整体结构，导致模拟结果精度较低。另一方面，三维立体数据虽然能够反映岩心的真实结构，但基于三维立体数据所构建的物理模型的数据存储量与运算量骤增，导致渗透率模拟周期延长、模拟效率和精度降低。

发明内容

有鉴于此，本发明的目的在于提供一种岩心渗透率模拟方法及装置，以解决现有技术的岩心渗透率模拟结果计算周期长、精度低的问题。

第一方面，本发明实施例提供了一种岩心渗透率模拟方法，包括以下步骤：

读取岩心的CT图像数据；

对CT图像数据进行二值化处理，生成岩心多孔介质的二值化CT图像数据；

对二值化CT图像数据进行三维重建，生成岩心的三维体数据；

通过三维体数据获取岩心的表征体元数据；

根据表征体元数据建立岩心的物理模型，并由物理模型计算宏观物理量，其中，物理模型分布函数的存储采用稀疏矩阵存储算法；

根据宏观物理量计算岩心渗透率。

物理模型分布函数采用稀疏矩阵存储算法进行存储的具体步骤为：

将流体格点的分布函数存储于形式为F[fluid_num][Q]的数组中，其中，每个流体格点相邻的Q-1个流体格点的索引位置存入辅助数组map[fluid_num][Q-1]中。

计算时，对于任一流体格点，先通过辅助数组获取该流体格点相邻格点的索引位置，然后根据索引位置将对应的分布函数写入F[fluid_num][Q]中，稀疏矩阵仅仅存储了流体格点的分布函数，避免了由于固体格点参与运算所导致的数据量、运算量都很庞大的问题，进而提高了物理模型的运算速度，特别是对于孔隙率较低的岩心，运算速度提高的更为明显。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，根据宏观物理量计算流经岩心中的介质的相对渗透率，即根据宏观物理量可以计算岩心渗透率和流经岩心的介质的相对渗透率。

结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，通过三维体数据获取岩心的表征体元数据，具体为：

在三维体数据中至少随机选取两个体素点，然后以每个体素点为中心选取边长为a的立方体，并统计每个立方体的孔隙度；

逐渐增加立方体的边长，使每个立方体的孔隙度均趋于稳定值的边长即为表征体元边长；

由表征体元边长所围成的立方体数据即为表征体元数据。

表征体元以最少的三维数据反映了岩心的整体结构，大大减少了参与模型计算的数据量，在减少运算量的同时提高了渗透率的计算精度，同时提高了计算效率。

结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，在根据表征体元数据建立岩心的物理模型之后，将介质在岩心中的接触角输入物理模型中，根据物理模型计算宏观物理量。

结合第一方面及其可能的实施方式，本发明实施例提供了第一方面的第四种可能的实施方式，其中，物理模型为多松弛动理学LBE模型，离散模型为D3Q19模型，多松弛动理学LBE模型的演化方程如下：

其中，I是单位张量，f_ai(x,t)是粒子密度分布函数，Ω_ai和分别是离散碰撞算子和离散的作用力项。

第二方面，本发明实施例还提供一种岩心渗透率模拟装置，包括：

读取模块，用于读取岩心的CT图像数据；

二值化模块，用于对CT图像数据进行二值化处理，生成岩心多孔介质的二值化CT图像数据；

三维重建模块，用于对二值化CT图像数据进行三维重建，生成岩心的三维体数据；

表征体元数据获取模块，通过三维体数据获取岩心的表征体元数据；

物理模型模块，用于根据表征体元数据建立岩心的物理模型，并根据物理模型计算宏观物理量，其中，物理模型分布函数的存储采用稀疏矩阵存储算法；

渗透率计算模块，用于根据宏观物理量计算岩心渗透率。

本实施例所提供的岩心渗透率模拟装置实现了岩心CT图像数据的读取、CT图像数据的二值化、三维重建，以及物理模型的创建和岩心渗透率的输出。

结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，通过渗透率计算模块计算流经所述岩心的介质的相对渗透率，即渗透率计算模块可以输出岩心的绝对渗透率和介质的相对渗透率，具有普适性，适用于不同介质相对渗透率的计算。

结合第二方面，本发明实施例提供了第二方面的第二种可能的实施方式，其中，表征体元数据获取模块具体用于：

在三维体数据中至少随机选取两体素点，然后以每个体素点为中心选取边长为a的立方体，并统计每个立方体的孔隙度；

逐渐增加立方体的边长，使每个立方体的孔隙度均趋于稳定值的边长即为表征体元边长；

由表征体元边长所围成的立方体数据即为表征体元数据。

本实施方式通过表征体元选取合适尺寸的数字岩心，既使岩石物理特性不再受岩石尺寸的影响，又将岩心数据量降到了最小，进而减少了渗透率计算的数据量，有助于提高渗透率的计算效率。

结合第二方面，本发明实施例提供了第二方面的第三种可能的实施方式，其中，还包括接触角模块，用于根据输入的接触角数值计算物理模型的宏观物理量。

结合第二方面及其可能的实施方式，本发明实施例提供了第二方面的第四种可能的实施方式，其中，物理模型为多松弛动理学LBE模型，离散速度模型为D3Q19模型，多松弛动理学LBE模型的演化方程如下：

其中，I是单位张量，f_ai(x,t)是粒子密度分布函数，Ω_ai和分别是离散碰撞算子和离散的作用力项。

本发明实施例带来了以下有益效果：

通过重建二值化的岩心CT图像数据生成岩心的三维体数据，通过三维体数据获取岩心的表征体元数据，使用表征体元数据建立物理模型，且采用稀疏矩阵存储算法存储物理模型的分布函数，然后将物理模型的输出值带入渗透率公式计算岩心绝对渗透率，表征体元数据在减少岩心三维体数据量的同时又能充分反应岩心的三维结构，进而能够输出准确的宏观量，在降低数据运算量的同时提高了岩心绝对渗透率的计算精度。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

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

图1为本发明实施例1提供的岩心渗透率模拟方法的流程图；

图2为本发明实施例1提供的表征单元体边长与孔隙度的关系图；

图3为本发明实施例1提供的不同空隙度时完全存储矩阵与稀疏存储矩阵的计算速度示意图；

图4为本发明实施例1提供的岩心绝对渗透率网格收敛性的示意图；

图5为本发明实施例1提供的理论接触角为150度时模型输出的预测接触角的三维示意图；

图6为本发明实施例1提供的理论接触角为150度时模型输出的预测接触角的二维示意图；

图7为本发明实施例1提供的理论接触角为90度时模型输出的预测接触角的三维示意图；

图8为本发明实施例1提供的理论接触角为90度时模型输出的预测接触角的二维示意图；

图9为本发明实施例1提供的理论接触角为45度时模型输出的预测接触角的三维示意图；

图10为本发明实施例1提供的理论接触角为45度时模型输出的预测接触角的二维示意图；

图11为本发明实施例1提供的岩心相对渗透率的数值模拟曲线；

图12为本发明实施例2提供的岩心渗透率模拟装置各模块连接示意图。

图标：

101-读取模块；102-二值化模块；103-三维重建模块；104-表征体元数据获取模块；105-物理模型模块；106-渗透率计算模块。

具体实施方式

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

目前岩心渗透率模拟结果的精度较低，基于此，本发明实施例提供的岩心渗透率模拟方法及装置可以提高岩心绝对渗透率以及介质的相对渗透率的精度。

为便于对本实施例进行理解，首先对本发明实施例所公开的一种岩心渗透率模拟方法进行详细介绍。

实施例1

本发明实施例提供的岩心渗透率模拟方法，如图1所示，包括以下步骤：

S101.读取岩心的CT图像数据

将真实岩心置于样品采集平台或是样品夹持器上，设置X射线强度、采集层厚、分辨率等参数，然后采集岩心的CT图像数据，本实施例中CT机分辨率优选为5.6938μm，然后读取所采集的二维CT数据。

S102.对CT图像数据进行二值化处理，生成岩心多孔介质的二值化CT图像数据

使用直方图阈值法或是迭代阈值法将采集的二维CT图像数据二值化，本实施例采用直方图法，且阈值优选阈值为146。

S103.对二值化CT图像数据进行三维重建，生成岩心的三维体数据

岩心的三维体数据信息量全，可以充分反映岩心的三维结构，克服了二维截面数据信息量少，带有片选偶然性的问题，从而使渗透率的计算结果具有很高的准确性和稳定性。

S104.通过三维体数据获取岩心的表征体元数据

表征体元(Representative Elementary Volume，简称REV)是指能够有效代表岩心宏观物理特性的最小岩心单元体。合适尺寸的数字岩心不仅使得岩石物理特性不再受岩石尺寸的影响，而且将岩心数据量降到了最小。本实施例中以岩石孔隙率变化作为REV尺寸的判断标准。

表征体元数据获取时，从三维体数据中至少随机选取两个体素点，如图2所示，本实施例中随机选择了三个体素点，然后以每个体素点为中心选取边长为a的立方体，并统计每个立方体的孔隙度；逐渐增加立方体的边长，使每个立方体的孔隙度均趋于稳定值的边长即为表征体元边长；由表征体元边长所围成的立方体数据即为表征体元数据。

需要说明的是，当立方体的边长增加时，孔隙度的变化值小于5×10^-4即认为孔隙度已趋于稳定，本实施例中阈值优选10^-4，且当边长a超过140后，孔隙度几乎趋于稳定值0.133，因此本实施例将140作为REV的基本尺寸，且此时岩心的孔隙度为0.133。

S105.根据表征体元数据建立岩心的物理模型，并由物理模型计算宏观物理量

物理模型为多松弛动理学LBE模型，离散速度模型为D3Q15或D3Q19模型，本实施例中的离散速度模型优选D3Q19，且多松弛动理学LBE模型的演化方程如下：

其中，I是单位张量，f_ai(x,t)是粒子密度分布函数，Ω_ai和分别是离散碰撞算子和离散的作用力项。

其中

F_a＝F_ap+F_ae＝-ρ_aK_aθ_a+ρ_ag_a

权系数为

离散速度为

分布函数对应的矩阵空间为

m_af_a＝(ρ_a,e,ε,j_x,q_x,j_y,q_y,j_z,q_z,3p_xx,3π_xx,p_ωω,π_ww,p_xy,p_yz,p_xz,m_x,m_y,m_z)

其中，ρ_a是组分a的质量密度；e和ε分别是能量和能量的平方；j＝j_x+j_t+j_z对应的是x，y和z方向上的动量；q＝(q_x+q_y+q_z)对应的是热流量；p_xx,p_xy,p_zx和p_ww是组分的压力张量；π_xx和π_ww对应的是粒子速度的四阶矩；m_x，m_y和m_z对应的是粒子速度的三阶矩。对角矩阵为

S＝diag(s₀,s₁,s₂,s₀,s₄,s₀,s₄,s₉,s₁₀,s₉,s₁₀,s₁₃,s₁₃,s₁₆,s₁₆,s₁₆,s₁₆)

其中，s₉＝s₁₃是粘性松弛时间，对应的粘性可以表示为

矩阵T的形式为

为了方便起见，在数值计算中将松弛时间都取为相同的值，并且不考虑流体内部的相互作用力。

另外，离散格子玻尔兹曼方法(Lattice Boltzmann Method，简称LBM)算法是一种显式时间推进算法，在每个时间步内，依次执行碰撞、迁移和边界处理三个步骤，多孔介质中有大量不参与计算的固体格点，本实施例中采用稀疏矩阵存储算法存储流体格点的分布函数，具体步骤为：将不参与计算的固体格点分离出来，而将流体格点的分布函数存储于形式为F[fluid_num][Q]的数组中，其中，每个流体格点相邻的Q-1个格点的索引位置存入辅助数组map[fluid_num][Q-1]中，本实施方式中流体格点的索引位置为流体格点的下标，在流体格点i执行k方向迁移操作前，先通过这个辅助数组查询其k方向的相邻格点在数组F中的位置map[i][k]，再执行写入操作，完成迁移过程，最后，在模拟计算结束并统计出宏观量后，再将其还原到原始的三维矩阵结构。

如图3所示，稀疏矩阵存储算法的使用大大减少了数据的存储量，尤其是孔隙率低的多孔介质，数据存储量大大减小，而且在计算过程中，也不需要遍历每个网格点进行判断，从而节省了计算时间，提高了数据的计算效率。

S106.根据宏观物理量计算岩心渗透率

绝对渗透率

采用步骤S105中的REV尺寸的数字岩心，其物理长度为L_x＝L_y＝L_z＝140×L_Res。其中，L_Res是CT扫描的最小分辨率，即5.6938um。数字岩心的孔隙率是0.133。采用常外力驱动，四周采用周期边界条件。外力加速度为g＝10^-4。其他物理量保持为dx＝dt，τ＝1.0。根据达西定律，离散格子玻尔兹曼方程(Lattice Boltzmann equation，简称LBE)绝对渗透率公式为

其中，Q是总体积流量，V是总体积，u是速度，ν是运动粘性，g是外力加速度。真实物理单位下的绝对渗透率K为

K＝K_LBE(L_Res)²

如图4所示，为了满足数值结果的网格无关性，本实施例分别模拟了网格数为N＝140，280，420，560的四种情况，当N＝560，岩心的绝对渗透率趋于平缓，可以排除网格数的影响，此时x、y和z方向的绝对渗透率的模拟值分别为K_x＝0.15D，K_y＝0.232D和K_z＝0.183D，且N＝560时，x方向的绝对渗透率的实验值为0.15137D，该实验值与模拟值吻合很好，有效地证明了LBE模型可以准确地预测多孔介质的单相流动。

需要说明的是，本实施例中的实验方法均采用以贝克莱-列维尔特前沿推进理论为出发点的非稳定方法，该方法相较于稳定法精度稍低，但稳定法对于致密岩心实验周期长，成本高，而非稳定法实验周期短，适用于各种岩心。

相对渗透率

多相流相对渗透率的计算涉及的影响因素较多，包括表面张力、粘性比以及润湿性等，本实施例所述模拟方法将可能影响计算精度的参数作为输入项输入到模型中。

本实施例模型还可以对接触角这一润湿特性的表征参数进行模拟，即接触角既作为输入项参与相对渗透率的计算，同时还可以根据输入的各个参数对接触角进行模拟，并将模拟结果作为输出项输出，具体参见附图5-10以及表1所示。

图5和图6示出了接触角输入值为150度时，本实施例模型输出的预测接触角的三维和二维示意图。

图7和图8示出了接触角输入值为90度时，本实施例模型输出的预测接触角的三维和二维示意图。

图9和图10示出了接触角输入值为45度时，本实施例模型输出的预测接触角的三维和二维示意图。

表1给出了理论接触角与模型输出的预测接触角的对比情况，从表中可以看出，45度接触角的误差为0.311％，90度接触角的误差为0.067％，150度接触角的误差为4.63％，均具有较高的准确率，因此，本实施例所述模型对于接触角的预测具有较高的稳定性。

表1理论接触角与预测接触角的对比表

本实施例中，相对渗透率的模拟以接触角为150度为例进行说明：

在相对渗透率测试实验中，岩心中注入水与试验油的粘度分别为1.27mpa·s和1.48mpa·s，粘度比接近1，通过弯钩法实验所得油滴与岩心壁面接触角为150度。为了计算方便，模拟时将粘性比设为1:1输入模型，同时将外力加速度为g＝10^-4、松弛时间为τ＝1.0、以及油滴与壁面接触角150度输入模型。

另外，为了保证网格无关性，网格数优选N＝560。油和水按照一定的饱和度S_w比例随机分布在多孔介质中，经过外力的驱动最终达到稳态。

模拟结果如图11所示，由油水两相的相对渗透率实验曲线可以看出，S_w＝0.3573时水的相对渗透率几乎为0。随着含水饱和度的增大，水相的相对渗透率不断增大，而油相的相对渗透率不断下降，在S_w＝0.55附近时，两个相对渗透率达到了相等的值，即等渗点。在数值模拟曲线中，S_w＝0.3573时，水的相对渗透率为0.016，随着含水饱和度的增大，水相的相对渗透率不断增大，而油相的相对渗透率不断下降，在S_w＝0.55附近时，两相对渗透率也达到了相等的值，即通过本实施例所述的多松弛动理学LBE模型可以模拟出具有较高精度的介质的相对渗透率。由此可见，表征体元数据的使用在保证岩心物理特性不变的情况下，减少了参与运算的原始三维数据，稀疏矩阵存储算法降低了运算过程中的数据存储与计算量，多松弛动理学LBE模型结合表征体元与稀疏矩阵存储算法可以在降低数据计算量的同时提高模拟结果的精确度和稳定性。

实施例2

如图12所示，本实施例提供了一种岩心渗透率模拟装置，包括：读取模块101，用于读取岩心的CT图像数据；二值化模块102，用于对CT图像数据进行二值化处理，生成岩心多孔介质的二值化CT图像数据；三维重建模块103，用于对二值化CT图像数据进行三维重建，生成岩心的三维体数据；表征体元数据获取模块104，通过三维体数据获取岩心的表征体元数据；物理模型模块105，用于根据表征体元数据建立岩心的物理模型，并根据物理模型计算宏观物理量；渗透率计算模块106，用于根据宏观物理量计算岩心渗透率。

该装置可以读取岩心的CT图像数据，然后对图像数据依次进行二值化、三维重建、表征体元数据获取以及建立物理模型，然后根据物理模型输出的宏观物理量输出岩心渗透率。

本实施例中的渗透率计算模块还能够计算流经岩心的介质的相对渗透率，本实施例用于计算油水的相对渗透率。

其中，表征体元数据获取模块104具体用于：在三维体数据中至少随机选取两体素点，本实施例中优选三个体素点，然后以每个体素点为中心选取边长为a的立方体，并统计每个立方体的孔隙度；逐渐增加立方体的边长，使每个立方体的孔隙度均趋于稳定值的边长即为表征体元边长；由表征体元边长所围成的立方体数据即为表征体元数据。

需要说明的是，当立方体的边长增加时，孔隙度的变化值在5×10^-4以内即认为孔隙度已趋于稳定，本实施例中的阈值范围为10^-4，如图2所示，当边长a超过140后，孔隙度几乎趋于稳定值0.133，因此本实施例将140作为REV的基本尺寸，且此时岩心的孔隙度为0.133。

物理模型为多松弛动理学LBE模型，离散速度模型为D3Q15或D3Q19模型，本实施例中的离散速度模型优选D3Q19，且多松弛动理学LBE模型的演化方程如下：

其中，I是单位张量，f_ai(x,t)是粒子密度分布函数，Ω_ai和分别是离散碰撞算子和离散的作用力项。

作为本实施例的优选实施方式，该装置通过稀疏矩阵存储模块存储流体格点的分布函数，具体为：采用形如F[fluid_num][Q]的稀疏矩阵存取流体格点的分布函数，其中，每个格点相邻的Q-1个格点的索引位置存入辅助数组map[fluid_num][Q-1]中，最后利用物理模型输出的宏观物理量计算最终的绝对渗透率和相对渗透率。

另外，本实施例中涉及的方法、公式如实施例1所述，在此不予赘述。

本发明实施例提供的岩心渗透率模拟装置，与上述实施例提供的岩心渗透率模拟方法具有相同的技术特征，所以也能解决相同的技术问题，达到相同的技术效果。

本发明实施例所提供的岩心渗透率模拟方法及装置的计算机程序产品，包括存储了程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

另外，在本发明实施例的描述中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (10)

1.一种岩心渗透率模拟方法，其特征在于，包括以下步骤：

读取所述岩心的CT图像数据；

对所述CT图像数据进行二值化处理，生成所述岩心多孔介质的二值化CT图像数据；

对所述二值化CT图像数据进行三维重建，生成所述岩心的三维体数据；

通过所述三维体数据获取所述岩心的表征体元数据；

根据所述表征体元数据建立所述岩心的物理模型，根据所述物理模型计算宏观物理量，其中，所述物理模型分布函数的存储采用稀疏矩阵存储算法；

根据所述宏观物理量计算所述岩心渗透率。

2.根据权利要求1所述的岩心渗透率模拟方法，其特征在于，根据所述宏观物理量计算流经所述岩心的介质的相对渗透率。

3.根据权利要求1所述的岩心渗透率模拟方法，其特征在于，通过所述三维体数据获取岩心的表征体元数据，具体为：

在所述三维体数据中至少随机选取两个体素点，然后以每个所述体素点为中心选取边长为a的立方体，并统计每个所述立方体的孔隙度；

逐渐增加所述立方体的边长，使每个所述立方体的孔隙度均趋于稳定值的边长即为所述表征体元边长；

由所述表征体元边长所围成的所述立方体数据即为所述表征体元数据。

4.根据权利要求1所述的岩心渗透率模拟方法，其特征在于，在所述根据所述表征体元数据建立所述岩心的物理模型之后，将所述介质在所述岩心中的接触角输入所述物理模型中，根据所述物理模型计算宏观物理量。

5.根据权利要求1-4任一项所述的岩心渗透率模拟方法，其特征在于，所述物理模型为多松弛动理学LBE模型，离散速度模型为D3Q19模型，所述多松弛动理学LBE模型的演化方程如下：

$f_{ai}(x+e_i\mathrm{\δ}t,t+\mathrm{\δ}t)-f_{ai}(x,t)={\mathrm{\Ω}}_{ai}+\mathrm{\δ}t{(I-\frac{1}{2}{T}^{-1}ST)}_{ij}{\stackrel{\‾}{F}}_{aj}$

其中，I是单位张量，f_ai(x,t)是粒子密度分布函数，Ω_ai和分别是离散碰撞算子和离散的作用力项。

6.一种岩心渗透率模拟装置，其特征在于，包括：

读取模块，用于读取所述岩心的CT图像数据；

二值化模块，用于对所述CT图像数据进行二值化处理，生成所述岩心多孔介质的二值化CT图像数据；

三维重建模块，用于对所述二值化CT图像数据进行三维重建，生成所述岩心的三维体数据；

表征体元数据获取模块，通过所述三维体数据获取所述岩心的表征体元数据；

物理模型模块，用于根据所述表征体元数据建立所述岩心的物理模型，并根据所述物理模型计算宏观物理量，其中，所述物理模型分布函数的存储采用稀疏矩阵存储算法；

渗透率计算模块，用于根据所述宏观物理量计算所述岩心渗透率。

7.根据权利要求6所述的岩心渗透率模拟装置，其特征在于，通过所述渗透率计算模块计算流经所述岩心的介质的相对渗透率。

8.根据权利要求6所述的岩心渗透率模拟装置，其特征在于，所述表征体元数据获取模块具体用于：

在所述三维体数据中至少随机选取两体素点，然后以每个所述体素点为中心选取边长为a的立方体，并统计每个所述立方体的孔隙度；

逐渐增加所述立方体的边长，使每个所述立方体的孔隙度均趋于稳定值的边长即为所述表征体元边长；

由所述表征体元边长所围成的立方体数据即为所述表征体元数据。

9.根据权利要求6所述的岩心渗透率模拟装置，其特征在于，还包括接触角模块，用于根据输入的接触角数值计算所述物理模型的宏观物理量。

10.根据权利要求6-9任一项所述的，其特征在于，所述物理模型为多松弛动理学LBE模型，离散速度模型为D3Q19模型，所述多松弛动理学LBE模型的演化方程如下：

$f_{ai}(x+e_i\mathrm{\δ}t,t+\mathrm{\δ}t)-f_{ai}(x,t)={\mathrm{\Ω}}_{ai}+\mathrm{\δ}t{(I-\frac{1}{2}{T}^{-1}ST)}_{ij}{\stackrel{\‾}{F}}_{aj}$

其中，I是单位张量，f_ai(x,t)是粒子密度分布函数，Ω_ai和分别是离散碰撞算子和离散的作用力项。