CN108267390B - 一种含纳米孔隙储层的气体渗透率确定方法 - Google Patents

Abstract

本发明涉及一种含纳米孔隙储层的气体渗透率确定方法，所述方法利用格子玻尔兹曼(LBM)模拟计算气体在岩的数字模型矩阵中的流动，统计得到不同进口压力条件下的平均流体速度

(m/s)；将不同进出口压力差所对应的平均流体速度用拟合目标函数进行最小二乘法拟合，从而得到岩未知的固有渗透率K_o(m²)和扩散系数D_k(m²/s)；然后求出滑脱因子b_k(Pa)，再求出岩的气体渗透率即表观渗透率K_a(m²)；该方法解决了无法测量过于致密的储层的渗透率这个难题。

Description

一种含纳米孔隙储层的气体渗透率确定方法

技术领域

本发明属于油气田开发技术领域，具体涉及一种含纳米孔隙储层的气体渗透率参数确定方法。先是通过CT扫描获得储层岩心内部结构状态，利用LBM数值模拟计算流量压差关系，进而确定储层岩心渗透率的方法。

背景技术

对于含纳米孔隙的储层来讲，孔道狭窄、渗透率极低是这种储层岩石的特有性质。这种储层岩石内孔隙直径范围在0.1nm-100nm之间，如此狭小的孔道，只有气体分子能够顺利流过，而一般液体是无法轻易流过。传统的柱塞实验由于仪器的精密程度有限，导致测量渗透率误差太大，无法确定页岩的气体渗透率。而且由于只有气体才能通过这种致密的孔隙，也不能使用压汞的方法来测量含纳米孔隙储层的孔隙直径和孔隙度，使得通过孔隙直径和孔隙度计算测量获得气体渗透率的方法也无法实现。

目前也有方法通过CT扫描含纳米空隙储层的内部结构，通过分析图片结果来大致估量储层的孔隙直径和渗透率，然后根据经验计算公式得到储层的气测渗透率。但是这种方法是根据原有的微米级孔隙的多孔介质所推算的，在含有纳米孔隙的储层中不再适用，所以这种方法所形成的误差也很大。

发明内容

本发明提供了一种含纳米孔隙储层的气体渗透率确定方法。该方法将含纳米孔隙的储层岩的CT扫描数据重构成数字模型，然后使用格子玻尔兹曼方法(Lattice BoltzmannMethod，简称为LBM)模拟类似柱塞实验的过程，根据测定的压力和流速关系来推算含纳米孔隙的储层岩石的气体渗透率。所述格子玻尔兹曼方法是一种数值模拟计算方法，可以用来模拟介观状态下纳米级多孔介质内气体流动状态,其优点在于可以大大减少实验过程中产生的误差，可以在模拟环境中实现实验时无法达到的条件，从而解决无法测量过于致密的储层的渗透率这个难题。

所述方法包括如下步骤：

(1)选定所述储层的岩石，并将其加工成合适尺寸及形状；

(2)将加工成型的岩石进行扫描，得到描述内部孔隙的数据；

(3)将数据进行处理，生成岩的数字模型矩阵，获取包括模型长度L(m)的数字模型矩阵参数；

(4)利用格子玻尔兹曼(LBM)模拟计算气体在岩的数字模型矩阵中的流动，确定环境压力P_a(Pa)，不断提高进口压力P_e(Pa)，然后根据格子玻尔兹曼(LBM)模拟结果，统计得到不同进口压力条件下的平均流速

(5)将不同进出口压力差所对应的平均流速

用拟合目标函数进行最小二乘法拟合，从而得到岩未知的固有渗透率K_o(m²)和扩散系数D_k(m²/s)；然后求出滑脱因子b_k(Pa)，再求出岩的气体渗透率即表观渗透率K_a(m²)；

其中，所述步骤(5)中将不同进出口压力差所对应的平均流速

用拟合目标函数进行最小二乘法拟合，从而得到岩未知的固有渗透率K_o(m²)和扩散系数D_k(m²/s)，具体通过以下公式得到：

其中μ是气体粘度，单位为Pa·s。

其中，所述扫描为CT扫描或电镜扫描。

其中，含纳米孔隙的储层为页岩储层。

其中，步骤(1)中加工成合适尺寸及形状具体为加工成5mm×5mm×5mm大小的立方形试件。

其中，滑脱因子根据以下公式得到：

其中μ是气体粘度，单位为Pa·s。

其中，步骤(5)中的表观渗透率K_a(m²)根据以下公式得到：

所述方法采用数值模拟和真实模型相结合的方法解决实验环境下无法测得页岩液测渗透率和扩散系数的问题。

附图说明

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

图1为LBM数值模拟方法的流程图。

具体实施方式

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

利用本发明可以应用于测量含纳米级孔隙的储层气体渗透率，而且模拟所测得实验数据符合常规认识。下面根据从现场取得的页岩岩的确定过程，结合附图进一步说明本发明。

由于页岩的平均孔隙直径在0.2-200nm范围内，其气测渗透率大约在10^-18m²这个量级。其致密的特性使得常规柱塞气测实验无法正常进行，即使在岩两端施加非常高的压力差也会因为压力无法稳定而不能有效的测量页岩的气测渗透率。所以本发明考虑到数值模拟可以使数据结果稳定可靠的优点下，将数值模拟LBM和CT扫描页岩内部孔隙结构结合起来，通过计算模拟来获得页岩岩的气测渗透率。其中CT扫描只是一种示例，还可以采用其他任何的扫描方式，比如电镜扫描等，具体步骤如下：

第一步：选定一块页岩岩石，并将其加工成合适尺寸及形状。为了真实反映页岩的气测渗透率，必须选用实际页岩岩石作为原始结构资料计算。根据CT扫描机的要求，将地下深层的页岩岩石加工成5mm×5mm×5mm大小的方形试件，这是为了方便CT扫描时完整成像，如果形状不规则，扫描出来的数据可能非常杂乱；如果模型太大，则精度会变得太低，无法观测到内部孔隙结构。

第二步：将加工成型的岩石进行CT扫描，尽量提高其扫描精度，一般页岩的平均孔隙直径在0.2-200nm范围内，所以必须将辨识精度提高到1μm以下，才能清晰的看到页岩内部比较大的孔道分布。比较小的孔道难以在CT扫描机中呈现，而且对页岩内气体渗流贡献甚小，所以可以忽略不计。这样可以得到关于这个岩的多帧扫描图像。

第三步：将数据代入计算机进行处理，生成岩的数字模型矩阵。所扫描的数据为0-1矩阵，而且只是每一断面的扫描数据，所以用MATLAB将扫描的数据还原为三维0-1矩阵，这样就可以得到一个有干湿相区分的数字模型矩阵，其中0代表湿相，1代表干相。湿相是指多孔介质中可以流过流体的空间，也是渗流时流体所占据的空间，而干相则指多孔介质中的固体部分。然后从中选取一个比较适合渗流通过的长方条矩阵A，作为下一步数值模拟需要的基本模型，记录这个模型的长宽高为L(m)、W(m)、H(m)，以及矩阵A的格子长度L'(无量纲)、W'(无量纲)、H'(无量纲)。用MTALAB统计模型内的孔隙格子数，即湿相的格子数V_pore(无量纲)，然后描线并记录所有孔隙的总长度L_pore(无量纲)，根据公式

可以计算出模型的平均孔隙直径

(无量纲)。这样所得的平均孔隙直径比用其他统计方法测得的平均孔隙直径会偏大，不过不会影响最终计算结果的精确度。

第四步：用开源LBM代码模拟计算气体在数字模拟岩中的流动。格子玻尔兹曼方法(LBM)是一种可以在介观下模拟纳微米级流体流动的数值方法。其优点在于计算量偏少，理论简单，操作简单易行，而且可以模拟介观离散流体。本专利中使用D3Q19模型进行数值模拟，D3Q19模型可以使用于模拟等温三维气体流动。参见图1，具体的LBM模拟过程如下：

步骤101：可以按照以下方法设置进行数值模拟的条件：

设置初始状态下真实密度ρ₀(kg/m³)为标准状态下的甲烷气体密度，而初始状态下的真实速度u₀(m/s)为0，即u₀＝0。格子初始密度矩阵ρ'₀(无量纲)为矩阵A的同长度三维矩阵，初始密度矩阵ρ'₀内所有元素都为1，即ρ'₀＝1；设置初始速度矩阵u'₀(无量纲)为矩阵A的同尺寸的三维矩阵，初始速度矩阵u'₀内所有元素都为0，即u'₀＝0。

真实物理密度矩阵ρ和格子密度矩阵ρ'、真实物理速度矩阵u和格子速度u'矩阵之间都存在着固定的比例关系，这个比例关系就是参考变量。参考变量是指实际物理变量和LBM中变量的比值，为了实现LBM变量与实际物理变量之间的转换，需要部分参考量：参考长度L_r(m)、参考密度ρ_r(kg/m³)、参考速度u_r(m/s)。

其中L、W、H、ρ和C_s为真实长、宽、高、密度和声速，而L'、W'、H'、ρ'、C_s'为LBM中的格子长、宽、高、格子密度和格子声速。实际物理量气体粘度μ(Pa·s)也是已知的，根据下述关系式：

在本模拟过程中，取真实长度L＝6μm、真实宽度W＝2.5μm、真实高度H＝2.5μm、真实环境压力P_a＝0.101MPa(即标准大气压值)、气体常数与温度乘积R_gT＝141933.671Pa·m³/kg，真实声速C_s＝376.743m/s、真实气体粘度μ＝1.1×10-⁵Pa·s；格子长度L'＝600、格子宽度W'＝250、格子高度H'＝250、初始状态下格子密度ρ'₀＝1、格子声速

所以初始状态下甲烷气体的真实密度为ρ₀＝P_a/R_gT＝0.7116kg/m³，而参考变量u_r＝652.538m/s、ρ_r＝ρ₀/ρ'₀＝0.7116kg/m³、L_r＝L/L'＝1×10-⁸m，格子粘度μ'＝μ/(L_rρ_ru_r)＝2.34。

步骤102：可以按照以下方法设置进行数值模拟的循环步骤：

然后设置进口处P_e＝0.5Mpa，则ρ_e＝P_e/R_gT＝3.5228kg/m³，ρ'_e＝ρ_e·ρ_r＝4.9505，这样格子初始密度矩阵ρ'₀在进口处截面的所有数值改为ρ_e，然后将新的密度矩阵ρ'进行迭代计算。其第一次迭代具体步骤为：

(1)根据下面的公式，将密度矩阵ρ'改为当前时间步内碰撞前的分布函数矩阵f(f为四维矩阵，无量纲)：

f＝ω_i·ρ'

其中，ω_i(无量纲)为D3Q19模型中的权函数。

(2)然后计算本迭代中的碰撞步，根据LBGK控制方程，即单松弛度离散玻尔兹曼控制方程，计算出当前时间步内经过碰撞后的分布函数矩阵。LBGK列于下方：

其中，f为碰撞前分布函数，f₁(无量纲)为碰撞后分布函数，f^eq(无量纲)为平衡态分布函数，τ(无量纲)为松弛时间。f^eq可以根据下列公式表达：

其中，u'为格子速度矩阵，e_i(无量纲)为微观速度分布。微观速度分布函数e_i和松弛时间τ配置如下：

这样就可以求得碰撞后分布函数f₁。

(3)然后进行迁移步，对f₁按照微观速度分布e_i进行平行迁移改变，得到迁移后的分布函数f₂(无量纲)，根据ρ'₁＝∑f₂，得到经过碰撞、迁移后的格子密度矩阵ρ'₁。

(4)再根据

可以得到经过碰撞、迁移后的格子速度矩阵u'₁。从格子速度矩阵u'₁中提取进出口的截面速度矩阵u'_e、u'_a(二维矩阵)，计算进出口速度的平均值u'_e、u'_a，并计算出进出口平均速度差值Δu'：

到此为止，第一次迭代完成。然后另下一时间步的碰撞前分布函数f替换为当前时间步内迁移后的分布函数f₂，并不断重复上述迭代，直到Δu'≈0，终止迭代。

上述为一次LBM计算模拟的完整过程，可以得到稳定状态下的出口平均流速u'_a，之后不断提高进口的压力边界条件，每次提高大约0.2MPa，一直到设置进口压力为2.1MPa结束，令平均流速

这样大概能得到一组共9对数据，9对数据只是示意性说明，当然也可以根据需要进行增加或减少，得到多对数据，能够实现拟合得出目标函数中的未知数即可。第五步：将不同进出口压力差P_e-P_a和所对应的平均流速

进行拟合。最小二乘法是常用的数值分析方法，简单有效，根据9对数据，可以拟合得出目标函数中的未知数，即固有渗透率Ko和扩散系数D_k，其中，将不同进出口压力差所对应的平均流速

用拟合目标函数进行最小二乘法拟合，从而得到岩未知的固有渗透率K_o(m²)和扩散系数D_k(m²/s)，具体通过以下公式得到：

其中，μ是气体粘度，单位为Pa·s。

然后根据滑脱因子的计算公式，可以得到克林伯格滑脱效应中使用的滑脱因子b_k。这个公式的优点是因为当渗流介质为纳微米级多孔介质时，传统的b_k计算方法和真实的数据误差太大，无法继续使用传统的b_k的计算公式，所以采用b_k＝μD_k/K_o这个计算公式。然后根据滑脱公式

可以求得页岩的气体渗透率，这样求得的页岩的气体渗透率不是一个定值，它是一个与进出口压力差有关的关系式。这比传统气测岩石渗透率更贴近理论。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种含纳米孔隙储层的气体渗透率确定方法，其特征在于，所述方法包括如下步骤：

(1)选定所述储层的岩石，并将其加工成合适尺寸及形状；

(2)将加工成型的岩石进行扫描，得到描述内部孔隙的数据；

(3)将数据进行处理，生成岩的数字模型矩阵，获取包括模型长度L，m的数字模型矩阵参数；

(4)利用格子玻尔兹曼LBM模拟计算气体在岩的数字模型矩阵中的流动，确定环境压力P_a，Pa，不断提高进口压力P_e，Pa，然后根据格子玻尔兹曼LBM模拟结果，统计得到不同进口压力条件下的平均流速

m/s；

(5)将不同进出口压力差所对应的平均流速

用拟合目标函数进行最小二乘法拟合，从而得到岩未知的固有渗透率K_o，m²和扩散系数D_k，m²/s；然后求出滑脱因子b_k，Pa，再求出岩的气体渗透率即表观渗透率K_a，m²；

其中，所述步骤(5)中将不同进出口压力差所对应的平均流速

用拟合目标函数进行最小二乘法拟合，从而得到岩未知的固有渗透率K_o，m²和扩散系数D_k，m²/s，具体通过以下公式得到：

其中，μ是气体粘度，单位为Pa·s；

其中，所述滑脱因子根据以下公式得到：

其中，μ是气体粘度，单位为Pa·s；

其中，所述步骤(5)中的表观渗透率K_a，m²根据以下公式得到：

2.根据权利要求1所述的方法，其特征在于，所述扫描为CT扫描或电镜扫描。

3.根据权利要求1或2所述的方法，其特征在于，含纳米孔隙的储层为页岩储层。

4.根据权利要求1-3之一所述的方法，其特征在于，步骤(1)中加工成合适尺寸及形状具体为加工成5mm×5mm×5mm大小的立方形试件。

Images