利用数字岩心抽象孔隙模型计算孔隙体积压缩系数

隋微波，权子涵，侯亚南，程浩然

（1.中国石油大学（北京）油气资源与探测国家重点实验室，北京 102249；2.中国石油大学（北京）石油工程学院，北京102249；3.深圳清华大学研究院，广东深圳 518057；4.清能艾科（深圳）能源技术有限公司，广东深圳 518057）

0 引言

岩石微观结构与其宏观物理性质密切相关。宏观性质或行为实质上是微观结构的集中体现，或者说微观结构控制宏观性质。Torquato[1]在研究非均质材料物理性质时表明，材料的宏观性质如弹性、渗透率等都会受到微观结构变化的影响。在细观力学研究领域，计算岩石弹性力学参数的理论方法一般有解析法和数值法两类。数值方法需要在利用场发射扫描电镜或CT设备观测岩石样品获得微观结构后，采用有限元等数值方法对岩石微观结构进行网格剖分，并模拟加载过程、计算应力应变场，再反推出弹性力学参数。比较具有代表性的工作是 Arns等[2]采用有限元方法计算Fountainebleau砂岩数字岩心的弹性参数，结果表明弹性参数的变化规律符合Gassmann理论。数值方法可以比较真实地重现岩石的微观孔隙结构，但是真正拟合孔隙结构时网格剖分难度很大，而且整体计算所需的时间成本非常高。解析方法一般基于Eshelby等效介质理论，将岩石中的孔隙抽象为球形、椭球形、圆柱形、扁裂缝与硬空心球壳等，采用解析方法计算岩石弹性力学有效性质。例如 Zimmerman[3]假设岩石中的孔隙为椭球形，根据不同的纵横比具体可以划分为球形孔、针状孔和币状孔，根据线弹性胡克定律以及 Eshelby等效介质理论等，推导出微观尺度的岩石变形理论。解析方法与数值方法相比可节约大量计算成本，但是由于一般不对岩心的具体微观孔隙结构进行刻画，大多用于理论方面的定性研究。

数字岩心技术是近年兴起的研究储集层微观孔隙结构、渗流特征参数和渗流机理的新手段[4-8]，因其对岩心样品尺寸要求低，实验及模拟计算可重复性强，成为非常规油气微纳米尺度条件下储集层性质和渗流机理研究的新途径[9-10]。虽然近年来数字岩心技术已形成了较为成熟的研究方法与理论体系，但是与真实的储集层情况和宏观油藏研究手段相比[11]，数字岩心技术在模拟微观渗流时一般没有考虑真实条件下储集层应力对微观孔隙结构和渗流情况的影响，忽略了应力敏感问题，这使得数字岩心技术在实际应用中还存在一定缺陷。本文研究数字岩心基本弹性参数（孔隙体积压缩系数）与微观孔隙结构的关系，为数字岩心的微观变形研究提供理论支持。在油藏工程中，孔隙体积压缩系数不仅是物质平衡方法计算中的重要输入参数[12]，同时也是渗流力学与试井分析中压力扩散方程的重要参数[13]。

本文采用Bentheimer砂岩等10种数字岩心，以数字岩心真实微观孔隙结构为基础，根据数字岩心微观孔隙结构分析孔隙形状参数分布特征，建立抽象孔隙模型，并采用解析方法结合等效介质理论，考虑微观结构变形的本构关系，求取数字岩心有效弹性模量和孔隙体积压缩系数，分析孔隙结构对孔隙体积压缩系数的影响。

1 岩石抽象孔隙模型的理论基础

1.1 单孔弹性变形

将孔隙抽象为椭球体最初是由 Sadowsky等[14-15]提出的，椭球体可以代表大部分孔隙形状，例如球体、针状圆柱体、薄片裂缝等。对于长半轴长度为A，短半轴长度为B的椭圆（见图1a），当椭圆沿长轴为中心进行空间旋转时得到长椭球体（Prolate Ellipsoid），也称为“针状孔”（见图 1b）；当椭圆以短轴为中心进行空间旋转时得到扁椭球体（Oblate Ellipsoid），也称为“币状孔”（见图1c）。

图1 长椭球、扁椭球纵横比定义示意图

椭球体纵横比（ρ）定义为不等长半轴与等长半轴之比，所以长椭球形的纵横比大于1（ρ＞1），扁椭球形的纵横比小于1（ρ＜1）。因此，通过求解球坐标系中的弹性问题，就能够获得其孔隙压缩性参数[3]。

对于考虑微观结构的多孔介质，Zimmerman认为对于微观尺度的孔隙来说，在多孔介质上所施加的外部应力相当于在无穷远处，因此可以只考虑正向应力，则孔隙体积应变可以表示为[3]：

若先不考虑孔隙空间中含有流体的情况，即孔隙压力不变，则上式变为

孔隙体积压缩系数是指应力发生变化后，孔隙体积的变化量与应力变化之前孔隙体积的比值，孔隙体积压缩系数Cpp可用来表征储集层岩石孔隙的变形尺度。下面首先分析单孔变形条件下孔隙体积压缩系数的解，分别考虑长、扁椭球孔变形和球形孔变形两种情况，然后将其引申为多孔变形的情况。

1.1.1 长、扁椭球孔变形

对于无限大、各向同性弹性介质中独立长椭球孔的受压变形问题，根据Zimmerman的单孔变形理论[3]，可以长椭球孔中心为原点建立直角坐标系，z轴为长半轴，再将直角坐标转换为长椭球坐标（见图2）：

其中l为长椭球长轴（z轴）上的焦距半长，下文推导中取l为单位长度。α取值范围为0到∞，β为0到π，γ为0到2π。当α为常数时，方程表示长椭球球体表面；β为常数时，方程表示双曲面；γ为常数时，方程表示过z轴的平面。球孔的表面相当于α=α0，两条等长的短轴为sinhα0，长轴为coshα0，纵横比ρ=cothα0。

图2 长椭球坐标转换示意图

为了方便表示，引入以下变量：

则该正交系下弧长表示为：

度量系数hi（i=α,β,γ）有以下关系：

由于所施加的压力是垂直于孔洞表面的单位压力，且孔洞无穷远处的应力为零，则边界条件为：

上述边界条件下的线弹性力学问题可以采用Sadowsky等[14]提出的 Bossinesq三函数法求解，从而获得位移（U）的解。通过进一步对位移正分量在整个孔洞表面的积分来计算孔隙体积的变化量，进而获得孔隙压缩系数。

假设孔隙所受应力变化为Δσ，则孔隙压缩系数表示为：

将（8）、（10）式代入（9）式，可以求得长椭球坐标系下的孔隙压缩系数为：

长椭球坐标系通过变换即可获得扁椭球坐标系，因此扁椭球的孔隙压缩系数可以直接由（11）式转换得到：

长椭球孔、扁椭球孔的体积应变可分别表示为：

1.1.2 球形孔变形

球形孔的体积应变为：

由（18）式可求得球形孔的孔隙体积压缩系数为：

1.2 多孔弹性变形

前面讨论了单个不同形状孔隙的压缩系数，表征为单个孔隙在单位压差下的体积变化量。整个孔隙体系的孔隙体积变化量则为每个独立孔隙体积变化量之和。但在真实岩石中，每个孔隙的应力场会因周围其他孔隙的存在而受到影响，周围孔隙的存在会使每个单独孔隙的孔隙压缩系数增加。

在真实多孔介质材料（如岩石）的弹性变形分析中，考虑到岩石的宏观力学性质（或传输性质）与其微观结构有着密切关系，常采用Eshelby[17]提出的等效介质理论，通过微观结构计算得到宏观材料的有效性质。本文将岩石等效为无限大基质中带有夹杂的混合物，即将孔隙看作其中的夹杂，骨架颗粒看作基质，则多孔介质岩石的有效弹性性质与骨架-孔隙的弹性模量、体积分数、微观几何形状、空间分布有关。求解多孔弹性变形条件下的孔隙体积压缩系数前，需要首先求解多孔条件下的岩石有效模量。

1.2.1 有效模量求解

多孔弹性变形条件下的有效模量可采用自洽法求取，该方法最早由Hill[18]和Budiansky[19]提出，是最常用的等效介质法。自洽法将复合材料中的夹杂（孔隙）视为球形，然后嵌入到未知有效模量的介质中进行有效性质的计算。Wu[20]将夹杂视为椭球形并提出了不同于球形夹杂的自洽模型。

本文考虑将任意形状和性质的夹杂（孔隙）都放入一个无限大的基质中，这个含有夹杂的基质其性质等于整个非均质材料的未知有效性质。假设夹杂和基质两相均为各向同性，同时考虑球形夹杂为任意分布，则有效体积模量和有效剪切模量计算公式为[21]：

该式为隐式方程，计算过程需要使用数值迭代法。对于非球形夹杂的情况，采用Wu[20]提出的自洽模型，将夹杂视为椭球形处理，则两相混合物模量的计算公式为：

1.2.2 孔隙体积压缩系数求解

考虑到真实储集层岩石均含有多种孔隙形状，可以根据不同纵横比分布的孔隙体积占比计算最终的孔隙体积压缩系数。设岩心样品孔隙纵横比分布区间为［ρ1,ρN］，其中纵横比为ρk的孔隙体积为Vk，孔隙体积占比为xk，则有：

混合孔隙压缩系数为：

2 数字岩心抽象孔隙模型建立

本文中选取的数字岩心样品涵盖了高、中、低渗透砂岩与碳酸盐岩，具体包括 Bentheimer砂岩、Doddington砂岩、Berea砂岩、Fountainebleau砂岩、Wilcox致密砂岩、Estaillades碳酸盐岩、大庆砂岩、塔里木砂岩、南海东部砂岩和新疆砂岩等 10种数字岩心样品。其中前 6种数字岩心样品来源于公开的数字岩心标准比对岩样数据体[22]，其他 4种分别是井下岩心（大庆和塔里木）和井下钻屑（南海东部和新疆）制样后CT扫描获得的数字岩心数据体。下面介绍利用Berea砂岩数字岩心建立抽象孔隙模型的方法，其他岩心的抽象孔隙模型建立方法与之相同。

本次研究中使用的 Berea砂岩数字岩心来源于英国帝国理工大学PERM研究课题组公布的网上数据[23]。该数据是将 Berea砂岩样品通过 CT扫描实验获得的1 024张1 024像素×1 024像素的二维灰度图像（分辨率为2.77 μm），使用三维成像软件将二维切片叠加得到岩样的三维灰度图像，同时采用分水岭算法进行阈值分割提取孔喉结构，获得二值化岩心图像。

基于数字岩心孔隙结构特征分析，应用 Avizo图像处理软件对孔隙整体区域进行分离，并对单个孔隙进行标注（见图 3，图中颜色用于区分孔隙个体）。应用惯性矩原理对所有标注孔隙进行参数分析，计算每个孔隙对应的具有相同标准的二阶中心矩椭圆的主轴长度、重心位置和欧拉角，通过Matlab编程将真实孔隙转换为具有相同特征参数的椭球体，则真实孔隙空间可转换为抽象孔隙模型（见图4，图中颜色用于区分孔隙个体）。

图3 Berea砂岩数字岩心抽象孔隙分离示意图

图4 Berea砂岩数字岩心孔隙抽象示意图

通过计算得到 Berea砂岩数字岩心样品抽象孔隙体积占比与椭球体纵横比之间的关系（见图 5），从该图中可以看出 Berea砂岩中各类孔隙纵横比分布范围为0.2～4.9。根据（27）式计算，可以由图中每个纵横比区间的孔隙体积占比计算最终的Berea砂岩孔隙体积压缩系数。其他9种数字岩心的孔隙体积压缩系数计算中也均需输入孔隙纵横比和孔隙体积占比的分布关系。

图5 Berea砂岩抽象孔隙体积占比与纵横比之间的关系

本文提出的数字岩心抽象孔隙模型方法适用于绝大部分砂岩和除鲕粒灰岩之外的大部分碳酸盐岩。主要原因是该方法的前提是假设岩心中的绝大部分孔隙为空间凸体，并在此基础上对椭球体进行抽象形成孔隙模型；而从Ketton灰岩的电镜扫描图片（见图6a）和 CT扫描获取的三维孔隙结构（见图 6b）可明显看出，鲕粒灰岩由于其成岩颗粒为圆形或椭圆形，颗粒间孔隙空间则呈明显凹体，因此该方法不适用。

图6 Ketton灰岩电镜扫描与CT成像获取的孔隙结构

3 孔隙体积压缩系数与孔隙结构

基于前文的岩石抽象孔隙模型理论基础和建立的数字岩心抽象孔隙模型，分别从单一类型孔隙和混合类型孔隙两方面研究孔隙体积压缩系数与孔隙结构的关系。

3.1 单一类型孔隙

假设岩心中只存在单一类型孔隙，考虑到岩石中发育多孔且孔隙纵横比不同，此处假设长椭球孔、扁椭球孔的纵横比符合高斯分布，且岩石骨架颗粒（石英）体积模量为37 GPa，剪切模量为44 GPa，岩石孔隙度为20%。在此基础上研究孔隙纵横比对孔隙压缩系数的影响。

首先考虑岩心中发育长椭球孔的情况，假设孔隙纵横比均值为 6.00～20.00，方差为 0.50～3.00，岩石孔隙度为 20%，则可得孔隙压缩系数Cpp,pro与纵横比均值、纵横比方差的关系（见图 7）。可以看出：①长椭球形孔的纵横比符合高斯分布时，孔隙压缩系数随纵横比均值的增大而增大，且纵横比小于10.00时增大速度较快；②方差越小孔隙压缩系数越大，说明在相同纵横比均值条件下，长椭球纵横比分布越集中孔隙压缩系数越大，在相同应力下的变形量也就会越大。

图7 长椭球形孔孔隙压缩系数与纵横比均值、方差的关系

同理，假设扁椭球的纵横比分布符合高斯分布，均值为0.10～0.80，方差为0.04～0.20，岩石孔隙度为20%，则可得孔隙压缩系数Cpp,ob与纵横比均值、纵横比方差的关系（见图8）。可以看出：①扁椭球形孔的纵横比按高斯分布时，孔隙压缩系数随纵横比均值的增大而减小，且方差越小减小速度越快；②与长椭球孔的结果相同，扁椭球孔的纵横比方差越小孔隙压缩系数越大，即扁椭球形孔纵横比分布越集中，孔隙压缩系数越大，在相同应力下的变形量就会越大；③当纵横比均值大于0.30时，孔隙压缩系数基本不随方差的变化而变化。

图8 扁椭球形孔孔隙压缩系数与纵横比均值、方差的关系

3.2 混合类型孔隙

真实岩心中极少存在单一类型孔隙，多以混合类型孔隙出现，其孔隙压缩系数可由（27）式计算。这里主要讨论不同类型孔隙体积占比、孔隙度对整体孔隙压缩系数的影响。

3.2.1 孔隙体积占比的影响

假设岩石孔隙度为 20%，长椭球形孔的纵横比为2.00，扁椭球形孔的纵横比为0.10，计算得混合孔隙压缩系数与不同类型孔隙（长椭球形、扁椭球形、球形）体积占比的关系（见图 9）。可以看出：①3种孔隙混合时，混合孔隙压缩系数与扁椭球孔隙体积占比呈正相关，与球形孔或长椭球形孔的孔隙体积占比呈负相关；②只有球形、长椭球形 2种孔隙混合时，孔隙压缩系数随球形孔的增多而增大。上述结论表明，在一定的纵横比前提下，长椭球形孔越多岩石越不容易变形，扁椭球形孔和球形孔越多岩石越容易发生变形，且扁椭球形孔的影响比球形孔更明显。

图9 不同混合模式下孔隙压缩系数与孔隙体积占比的关系

3.2.2 孔隙度的影响

假设长椭球形孔的纵横比为 2.00，扁椭球形孔的纵横比为 0.50，按 4种模式进行混合：①长椭球孔与扁椭球孔的混合比为1∶1；②扁椭球孔与球形孔的混合比为1∶1；③长椭球孔与球形孔的混合比为1∶1；④长椭球孔、扁椭球孔、球形孔的混合比为1∶1∶1。根据（27）式可以获得孔隙压缩系数与孔隙度的关系（见图10）。可以看出，孔隙压缩系数随孔隙度的增加而增大，但不同的混合模式下压缩系数及其增长速度有较大差异：①混有扁椭球孔的 3种模式，相同孔隙度条件下，孔隙压缩系数远大于长椭球孔与球形孔的混合模式，且孔隙度越大相差幅度越明显，说明具有扁椭球孔的岩石更易发生变形，混合孔隙压缩系数更大；②长椭球孔、扁椭球孔、球形孔同时存在时孔隙压缩系数小于长椭球与扁椭球孔混合、球形孔与扁椭球孔混合的情况，说明 3种孔隙同时存在会增加岩石的抗压能力；③扁椭球孔与球形孔混合，其孔隙压缩系数比扁椭球孔、长椭球孔混合的孔隙压缩系数略大，说明球形孔比长椭球孔更容易发生变形。

图10 不同混合模式下孔隙压缩系数与孔隙度的关系

3.3 真实数字岩心孔隙

图11 9种真实数字岩心样品孔隙分布图

表1 10种数字岩心抽象孔隙模型中不同类型孔隙纵横比与体积占比情况

除Berea砂岩外的 9种真实数字岩心样品孔隙分布如图 11所示（图中颜色用于区分孔隙个体）。在前述理论方法的基础上，对这10种真实数字岩心建立了抽象孔隙模型并计算了不同类型孔隙体积占比和纵横比分布范围（见表1）。可以看出，10种典型数字岩心中长椭球孔、扁椭球孔纵横比分布范围分别为 1.11～12.75和0.11～0.90。由于岩心最终的孔隙体积压缩系数由不同类型孔隙的体积占比计算，因此不同类型孔隙的体积占比与数量占比相比更具有比较意义。根据前文理论分析结果，对于岩心孔隙压缩系数影响较大的几个因素包括各类孔隙体积占比、纵横比分布范围、均值和方差。

根据前述理论模型、表 1中孔隙结构特征参数计算10种数字岩心样品的孔隙体积压缩系数（见图12），由图中可看出，孔隙体积压缩系数总的变化趋势是随孔隙度增高而增大。其中孔隙度较大的3个样品中（表中序号第1—3）Bentheimer和Doddington砂岩样品孔隙特征参数整体很接近，孔隙体积压缩系数受宏观孔隙度差异影响有所不同，Doddington砂岩孔隙度稍大因此孔隙体积压缩系数较大；Estaillades样品与Doddington样品孔隙度相同，但是 Estaillades样品具有更多的球形孔和长椭球孔，因此压缩系数相对较小。孔隙度中等的5个样品（表中序号第4、7—10）中，Berea砂岩样品比其他4个具有更多的球形孔，长椭球孔隙纵横比均值较小且分布范围非常小，扁椭球孔隙纵横比均值较大，因此压缩系数相对较低。孔隙度最低的 2个样品虽然孔隙压缩系数相差不大，但是其变化规律与其他样品略有不同，Wicox砂岩中长椭球孔体积占比高于Fountainbleau砂岩，扁椭球孔体积占比低于Fountainbleau砂岩，如果只从不同孔隙占比的角度分析，Fountainbleau砂岩应该具有更大的孔隙压缩系数，但是这与模型的计算结果不同。通过对数据的深入分析发现，这一结果与两种特低孔隙度样品中扁椭球孔纵横比的最小值有关，Wicox砂岩虽然扁椭球孔含量相对较低，但是样品中纵横比在 0.1～0.2的扁椭球孔更多，Fountainbleau砂岩中扁椭球孔纵横比最小值是0.22，这一区别最终造成了与Fountainbleau砂岩相比稍高的孔隙压缩系数。这一影响因素在以后的研究中将进行更加系统的分析。

图12 10种真实数字岩心孔隙压缩系数与孔隙度的关系

4 结论

数字岩心的孔隙体积压缩系数与孔隙度、孔隙纵横比以及不同类型孔隙体积占比有关：①长椭球孔孔隙压缩系数与孔隙纵横比正相关，扁椭球孔孔隙压缩系数与孔隙纵横比负相关。②孔隙纵横比满足高斯分布且纵横比均值相同，长椭球孔与扁椭球孔纵横比分布越集中孔隙压缩系数越大，相同应力条件下变形量越大。③孔隙压缩系数随孔隙度增大而增大。④当岩心中发育多种孔隙类型时，孔隙压缩系数与不同类型孔隙体积占比相关，孔隙度一定时，扁椭球孔、球形孔越多岩石越容易发生变形，孔隙压缩系数越大；长椭球孔越多岩石越不容易变形，孔隙压缩系数越小。

采用孔隙体积压缩系数解析计算方法计算10种典型数字岩心样品孔隙体积压缩系数，所得结果与理论模型所述规律相符，该方法可用于数字岩心样品的孔隙体积压缩系数计算。

符号注释：

A——椭圆长半轴长度，m；Ai——求解位移积分时的常数项，无因次；B——椭圆短半轴长度，m；Cpp——单孔孔隙体积压缩系数，Pa−1；——混合孔隙体积压缩系数，Pa−1；Cpp,ob——扁椭球孔隙体积压缩系数，Pa−1；Cpp,pro——长椭球孔隙体积压缩系数, Pa−1；Cpp,s——球形孔隙体积压缩系数，Pa−1；——纵横比为kρ的孔隙体积压缩系数，Pa−1；——含扁椭球孔隙的多孔介质孔隙体积压缩系数，Pa−1；——含长椭球孔隙的多孔介质孔隙体积压缩系数，Pa−1；G——剪切模量，Pa；G0——骨架相剪切模量，Pa；G1——孔隙相剪切模量，Pa；——有效剪切模量，Pa；h，hα，hβ，hγ——度量系数，m−1；i，j——椭球坐标系坐标轴编号；K0——骨架相体积模量，Pa；K1——孔隙相体积模量，Pa；K*SC——有效体积模量，Pa；l——长椭球长轴（z轴）上的焦距半长，m；N——孔隙纵横比分布的区间数目；p——压力，Pa；P，——椭球坐标系中间变量，m；Q，——椭球坐标系中间变量，m；Q0——椭球坐标系长半轴长度，m；——椭球坐标系短半轴长度，m；r——球形孔半径，m；Rpro，Rob——中间变量，无因次；s——椭球坐标系弧长，m；SG——计算有效剪切模量的孔隙几何因子，无因次；SK——计算有效体积模量的孔隙几何因子，无因次；U——位移，m；Vk——纵横比为kρ的孔隙体积，m3；Vp——孔隙体积，m3；xk——纵横比为kρ的孔隙体积占比，无因次；x，y，z——直角坐标系坐标，m；α，β，γ——椭球坐标系坐标，（°）；0α——椭球坐标系中椭球孔表面处的α坐标值，（°）；Bα——Biot系数，无因次；δm,n——克罗内克函数，m=n时为1，m≠n时为0，无因次；εp——孔隙体积应变，无因次；εp,pro，εp,ob，εp,s——长椭球、扁椭球、球形孔孔隙体积应变，无因次；ξ——孔隙表面积，m2；ρ——纵横比，无因次；1ρ——纵横比最小值，无因次；kρ——纵横比分布区间内特定值，无因次；Nρ——纵横比最大值，无因次；ρob——扁椭球孔纵横比，无因次；ρpro——长椭球孔纵横比，无因次；σ——正向应力，Pa；Δσ——应力变化，Pa；υ——泊松比，无因次；φ——孔隙度，%。