孔隙尺度各向异性与孔隙分布非均质性对多孔介质渗透率的影响机理

李滔，李闽，荆雪琪，肖文联，崔庆武

（1.西南石油大学油气藏地质及开发工程国家重点实验室，成都 610500；2.中石化新星（北京）新能源开发有限公司四川分公司，成都 610500；3.中石化中原石油工程有限公司钻井一公司，河南濮阳 457000）

0 引言

多孔介质的固有渗透率（文中简称渗透率）由其孔隙结构决定[1-2]。在油气开采领域，储集层渗透率对储集层评价和开发有着重要影响[3]。学者采用解析法[4-5]、实验法[6-10]和模拟法[11-14]，建立了多种渗透率与孔隙结构参数的关系式，多孔介质渗透率主要与孔隙度（或连通孔隙度）、孔径、迂曲度、比表面积和固相颗粒大小等有关。现有渗透率模型存在较大差异，且只适用于各向同性、相对均质的常规储集层[15]。与常规储集层相比，非常规储集层（致密砂岩、页岩等）岩石的孔隙结构具有以下特点：①微纳米级孔隙发育[16]，孔隙形态复杂[17-20]；②孔隙尺度各向异性强[21-23]；③孔隙分布非均质性强[17,24-25]。孔隙尺度各向异性与孔隙分布非均质性均显著影响多孔介质渗透率[22-23,26-27]，是非常规储集层不同于常规储集层的重要特征。Gu等[28]通过中子散射提出孔隙尺度各向异性的计算模型，但模型中包含了需测量的中子散射强度，极大限制了该模型的可应用性。随后，Wang等[26]基于形态学原理提出了具有较强适用性的多孔介质各向异性和孔隙分布非均质性计算模型。

目前对多孔介质各向异性和孔隙分布非均质性的定量研究相对缺乏，孔隙尺度各向异性与孔隙分布非均质性影响多孔介质渗透率的机理认识不清[21]，基于微CT图像重构三维数字岩心定量评价储集层孔隙尺度各向异性和孔隙分布非均质性的研究也未见文献报道。随着成像技术的发展，利用物理实验方法可直接观测多孔介质的微纳米级孔隙及其连通关系。采用四参数随机生成（QSGS）算法能够较高效地构建孔隙尺度各向异性、孔隙分布非均质性强的多孔介质，且能反映天然多孔介质绝大部分的形态复杂性[29]。同时，采用格子-玻尔兹曼方法（LBM），可直观、高效地处理流体内部及流-固间的相互作用[30]，实现多孔介质孔隙尺度的模拟[26,31]。这些均为非常规储集层渗透率影响机理的研究提供了条件。本文结合微CT扫描实验、QSGS算法和LBM，进一步揭示孔隙尺度各向异性与孔隙分布非均质性影响多孔介质渗透率的机理。

1 致密砂岩微观孔隙结构性质

1.1 三维数字岩心的建立

①四川盆地某气藏的4块致密砂岩岩心（H6-1、H6-2、H106-1和H106-2）的铸体薄片、扫描电镜和X-衍射实验结果表明，其孔隙空间由晶间微孔隙、粒间孔隙、粒内/粒间溶蚀微孔隙和粒间微缝等组成；成分以岩屑长石砂岩、长石岩屑砂岩为主，多为中粒、细—中粒砂岩，次为中—粗砂岩，少量为粉砂岩；胶结物主要为伊蒙混层、伊利石和黏土杂基。

②从4块岩样的端面钻取近似小圆柱（实际物理尺寸见表1），并对其进行微CT扫描（扫描测试电压150 kV，电流90 μA，扫描时间120 min，图像分辨率620 nm），扫描次序与岩样轴向一致。采用专业软件对岩样二维CT图像进行中值过滤和二值化处理，准确获取岩石骨架和孔隙空间，随后经图像切割完成微观结构成像。图1展示了H6-1和H106-1岩样单元体（800像素×800像素×800像素，对应实际物理尺寸0.496 mm×0.496 mm×0.496 mm）二值化前后以及孔隙空间的三维图像。可以看到致密砂岩孔隙形态复杂且极不规则。

表1 柱体岩样实际物理尺寸

③对单元体进行数值重构，建立致密砂岩的三维数字岩心，并计算得到4个单元体的孔隙度分别为2.60%、4.49%、4.51%和2.61%。

1.2 各向异性和非均质性的定量评价

基于已建立的三维数字岩心，采用Wang等[26]提出的多孔介质各向异性模型：

和孔隙分布非均质性模型：

定量评价致密砂岩的孔隙尺度各向异性与孔隙分布非均质性。

孔隙尺度各向异性因子被定义为多孔介质中所有孔隙的宽度之和与所有孔隙高度之和的比值，且各向异性强弱与A值成正比（A=1时表示各向同性）；孔隙分布非均质性因子表示单元块孔隙度的相对标准偏差（取决于观测尺度，即单元块大小），非均质性的强弱与H值成正比。

图1 基于CT图像获取的致密砂岩岩心三维图像

图2 孔隙尺度各向异性与数字岩心立方体边长的关系

图3 孔隙分布非均质性与数字岩心立方体边长的关系

选取线条间隔为5个格子，单元块长、宽、高均为25个格子，计算得到4块致密砂岩的孔隙尺度各向异性和孔隙分布非均质性与数字岩心立方体边长的关系（见图2、图3）。图中可以看出尽管4块致密砂岩三维数字岩心的孔隙度存在一定差异，但孔隙尺度各向异性、孔隙分布非均质性的变化却较为一致。图2显示，当立方体边长达到300个像素点时，4块致密砂岩的孔隙尺度各向异性趋于稳定，x、y和z方向（以岩样轴向为x，水平面内垂直于x的方向为y，垂直于x、y构成平面的方向为z，建立空间直角坐标系）的A值分别为1.38～1.55、0.70～0.77和0.91～1.02，岩样轴向的各向异性因子最大（这里将各向异性因子最大的方向定义为主流方向），其余2个方向相差不大。图3显示，立方体边长达到600个像素点，致密砂岩孔隙分布非均质性趋于稳定，H值范围为3.29～3.74。

2 三维随机多孔介质的LBM模拟

由以上分析可知，微CT图像能够较准确地反映岩样孔隙空间的三维形貌以及孔隙间的配置情况，但受分辨率限制，微CT扫描通常难以完整获取非常规储集层中的纳米级孔隙，由此建立的数字岩心孔隙空间不连通。对此，可在致密砂岩孔隙尺度各向异性、孔隙分布非均质性分析结果的基础上，采用QSGS方法构建三维各向异性、非均质性多孔介质模型，并运用三维19个离散速度方向（D3Q19）的多弛豫时间格子-玻尔兹曼（MRT-LB）模型模拟其中流体的运移，研究对渗透率的影响机理。

2.1 多孔介质构建

2.1.1 QSGS方法

QSGS方法构建三维多孔介质的具体过程[29]可分为3步：①设定构造区域与随机分布固相生长核（Cs），Cs的概率（Pcd）不能大于最终构建三维多孔介质的孔隙度（φ），即Pcd≤φ；利用随机数发生器生成[0，1]间的随机数，随机数小于Pcd的节点被指定为固相生长核；②给定不同方向的生长速度（Pk，其中k代表固相生长方向），利用随机数发生器生成[0，1]间的随机数，当Cs节点的随机数小于Pk时，则沿k方向的邻点生长；共设计26个生长方向，可分为侧线方向（6个）、对角线方向（12个）和直径线方向（8个）；当各方向的生长速度相同时，可获得各向同性多孔介质；当各方向的生长速度不同时，可获得各向异性多孔介质；③重复步骤②直到达到设定的多孔介质孔隙度。

为更好地研究孔隙分布非均质性对三维多孔介质渗透性的影响，可采用2个相互独立的生成过程来构建多孔介质[26]：①细结构生成较小的固相，固相生长核概率为Pr，孔隙度为Rφ；②粗结构生成较大的固相，固相生长核概率为Pc，且Pc远小于Pr；2个生成过程相结合，任一过程占据的节点均被设为固相，直至达到设定的孔隙度。

2.1.2 孤立孔隙删除

QSGS方法构建的三维多孔介质中包括孤立孔隙（或称死孔隙）和连通孔隙，连通孔隙是流体运移的通道。对四连通区域标记算法[32]进行改进可以去除三维多孔介质中的孤立孔隙：①扫描三维多孔介质中的孔隙和固相，将孔隙标记为0；②检查每个孔隙像素点东、南、西、北、上和下6个方向的邻近孔隙像素点，如果邻近的孔隙像素点均未被标记，则给该像素点设定一个新的标记数，新标记数从1开始递增，且不重复；如果只有一个邻近孔隙像素点被标记，则将孔隙像素点设定为与其相同的标记数；如果两个或以上的邻近孔隙像素点已被标记，则将孔隙像素点设定为已标记邻近孔隙像素点中最小的标记数；③检查三维多孔介质中所有的孔隙像素点，将相互连通的相邻孔隙像素点均使用其中最小的标记数标记；④检查各标记数所代表的孔隙空间是否同时与三维多孔介质出入口端相连。

2.2 LBM模型及验证

与单松弛（SRT）模型相比，多松弛（MRT）模型的碰撞算子中引入了多个松弛时间，计算精度更高且稳定性更好，也更适用于复杂多孔介质渗流问题的研究。MRT-LB模型的演化方程为[30]：

式中M是19×19的转换矩阵[33]，M-1是M的逆矩阵，可由高斯列主元消去法求得，S为对角矩阵：

式中非零矩阵元表示不同物理量（如：能量矩、能量平方矩、能量通量矩和应力张量矩等）的弛豫速率，取值分别为[34]：

D3Q19模型中，平衡态分布函数可表示为[34]：

式中宏观流体的密度、速度和运动黏度可分别表示为：

模拟中流体流动均为层流（Re值较小），结合达西定律可推导得三维多孔介质的渗透率计算公式[1]：

三维多孔介质的迂曲度等价计算公式[35]为：

为验证以上LBM算法的准确性，运用D3Q19 MRT和D3Q19 SRT模型（τR取1）分别模拟了体心立方中的流体运移。体心立方的孔隙度设为58.2%，边长设为320 nm，网格数分别为16×16×16、32×32×32、48×48×48、64×64×64、80×80×80、96×96×96和100×100×100。采用（8）式计算渗透率，并将结果与体心立方的渗透率解析解[36]作对比（见图4）可知，MRT模型模拟结果的准确性高于SRT模型；由于体心立方中球体半径与理论值始终存在一定的偏差[26]，随着网格数的增加，模拟结果只能越来越接近解析解而不能完全一致，即使在16×16×16网格条件下，MRT模型模拟结果与解析解的相对偏差最大仅为4.09%，说明MRT模型可准确模拟三维多孔介质中的流体运移。

图4 LBM模拟渗透率与解析解结果对比

2.3 模拟结果

2.3.1 三维孔隙尺度各向异性多孔介质的模拟

结合致密砂岩孔隙尺度各向异性的分析结果，同时为避免孔隙形态显著变化对模拟结果的影响，将x与y、z方向设定为不同的固相生长速度（其中y和z方向相等），采用QSGS算法构建三维孔隙尺度各向异性多孔介质（基础参数见表2）。为消除孔隙结构随机性对多孔介质渗透率的影响，同一组参数均构建10个随机多孔介质模型，共计180个。

表2 三维各向异性多孔介质模型QSGS算法与MRT-LB模拟中的基础参数

删除其中孤立孔隙，将线条间隔设为4个格子，运用（1）式计算多孔介质的孔隙尺度各向异性因子，得到多孔介质模型连通孔隙度、比表面积与孔隙尺度各向异性的关系（见图5、图6）。由图可知，多孔介质的各向异性因子介于0.98～1.62；各向异性因子由0.98增加至1.62时，多孔介质比表面积增加0.33%～2.55%，连通孔隙度减小1.31%～2.00%；随着Pcd的增加，多孔介质中小孔隙增多，比表面积显著增加（3种Pcd取值对应的比表面积均值分别为4.1×105cm2/cm3，4.7×105cm2/cm3，5.3×105cm2/cm3）；连通孔隙度受Pcd的影响相对较小，与Pcd呈负相关。

图5 各向异性模型比表面积与各向异性因子的关系

图6 各向异性模型连通孔隙度与各向异性因子的关系

采用D3Q19 MRT-LB模型模拟三维各向异性多孔介质中流体运移（基础参数见表2），模拟过程中的模拟偏差可表示为：

当模拟偏差小于10-6时，认为模拟结果不再变化，结束迭代过程。

图7为Pcd=0.002 5，A=1.55时各向异性多孔介质的模拟结果。以各向异性多孔介质的模拟结果为基础，采用（8）式、（9）式分别计算三维各向异性多孔介质模型的渗透率和迂曲度。分析渗透率、各向异性、迂曲度间的相互关系（见图8—图10）可知，各向异性多孔介质模型的渗透率与各向异性因子呈正相关；渗透率与迂曲度负相关；迂曲度与孔隙尺度各向异性因子呈负相关；各向异性因子由0.98增加至1.62时，Pcd为0.002 5，0.005 0和0.010 0的多孔介质渗透率分别增加80.20%，60.97%和57.38%，迂曲度分别减少17.72%，11.02%和6.06%。

图7 三维各向异性多孔介质模拟结果

图8 各向异性模型渗透率与各向异性因子的关系

图9 各向异性模型渗透率与迂曲度的关系

图10 各向异性模型迂曲度与各向异性因子的关系

这里需要说明的是：①数字岩心模拟计算多孔介质渗透率所得结果与单位网格长度密切相关，单位网格长度增加1个数量级时，多孔介质的渗透率至少增加2个数量级；②非常规储集层纳米级孔隙发育，为更好研究孔隙的微观结构及其与渗透率的相关性，模拟中选取了较小的单位网格长度（10 nm），渗透率模拟计算结果与实际岩心的测量值差异较大，但渗透率与孔隙微观结构的相关性与单位网格长度的取值无关。

2.3.2 三维孔隙分布非均质性多孔介质的模拟

采用QSGS算法生成三维非均质多孔介质模型（基础参数见表3）。这里x、y和z方向均设定相同的固相生长速度，忽略各向异性的影响，同时网格数和单位网格长度保持与三维各向异性多孔介质模型一致。为消除孔隙结构随机性的影响，同一组参数下均构建10个随机多孔介质模型，共计100个。构建完成后删除其中孤立孔隙，将单元块的长、宽和高均设为16个格子，采用（2）式计算多孔介质孔隙分布非均质性因子，分析连通孔隙度、比表面积与孔隙分布非均质性的关系（见图11、图12）可知，不同细结构孔隙度多孔介质的孔隙分布非均质性因子的计算结果分布范围有差距，当Rφ=40%时，孔隙分布非均质性因子相对较小，计算结果分布范围为0.89～1.13；当Rφ=60%时，非均质性较强，计算结果分布范围为1.03～1.42；多孔介质的连通孔隙度受非均质性的影响较弱，比表面积与非均质性呈现明显的负相关关系；当Rφ分别取40%和60%时，非均质性因子在其分布范围内增加，多孔介质连通孔隙度变化幅度均小于4.57%，而比表面积变化幅度较大，分别减小了21.55%和17.76%。

表3 三维非均质多孔介质模型QSGS算法中的基础参数

图11 非均质模型连通孔隙度与非均质性因子的关系

图12 非均质模型比表面积与非均质性因子的关系

运用D3Q19 MRT-LB模型模拟非均质多孔介质中流体运移（基础参数见表2）。以非均质模型的模拟结果为基础，运用（8）式和（9）式分别计算三维非均质性多孔介质的渗透率和迂曲度，分析渗透率、非均质性、比表面积、迂曲度间的相互关系（见图13—图16）可知，渗透率与孔隙分布非均质性因子、细结构孔隙度呈正相关[26]；Rφ越大，多孔介质中较大的孔隙越多，渗透率越大；当Rφ分别取40%和60%时，非均质性因子在其分布范围内增加，多孔介质渗透率分别增加136.01%和344.53%；渗透率与比表面积无明显相关性；渗透率与迂曲度呈明显的负相关；迂曲度与非均质性因子无明显的相关性。

3 分析与讨论

图13 非均质模型渗透率与非均质性因子的关系

图14 非均质性模型渗透率与比表面积的关系

图15 非均质性模型渗透率与迂曲度的关系

三维随机多孔介质模拟中，三维各向异性多孔介质中固相颗粒量级均为单位网格大小[37]，而非均质性多孔介质中，粗结构固相生长核概率减小，固相颗粒的大小逐渐增加。在此基础上，主要研究了渗透率、比表面积、迂曲度、连通孔隙度等物性参数与孔隙尺度各向异性、孔隙分布非均质性的关系。

这里进一步研究孔隙尺度各向异性、孔隙分布非均质性对随机多孔介质孔隙大小和分布的影响，以及各向异性、非均质性对复杂多孔介质渗透率的影响机理。采用最大球法获取多孔介质的孔隙直径、不同孔隙占孔隙空间的体积比例。最大球法获取三维多孔介质孔隙特征参数的步骤包括：①确定最大球半径的相关概率；②寻找最大球；③删除冗余球（包括单个及多个包含关系）等。即使多孔介质中孔径大小和分布均相同，孔隙间的配置关系也会导致渗透率的改变[1]，为消除这一影响，同一组基础参数均构建了10个随机多孔介质模型。

图16 非均质性模型迂曲度与非均质性因子的关系

多孔介质中流体的运移主要发生在较大的连通孔隙中，因此采用最大球法抽取时不考虑只包含单个节点的孔隙。图17、图18为各向异性（Pcd=0.002 5）、非均质性（Rφ=60%）多孔介质模型的孔径分布，横坐标为孔隙直径（以网格数量表示），纵坐标为同一组基础参数下10个随机多孔介质模型孔径分布中对应孔隙体积所占比重的均值。分析图17可知，不同生长速度（对应各向异性）条件下多孔介质的孔径分布形态十分接近，相对偏差为1.78%～8.64%；图18显示出不同粗结构固相生长核概率（对应非均质性）条件下多孔介质的孔径分布形态也十分接近，相对偏差为1.82%～6.09%；相对偏差大小与各向异性或非均质性强弱无关。

图17 各向异性模型的孔径分布

图18 非均质性模型的孔径分布

多孔介质的孔隙结构参数与孔隙尺度各向异性、孔隙分布非均质性的关系（见表4）可归纳为：①三维各向异性多孔介质模型中，Pcd相同，其连通孔隙度、孔径、比表面积和固相颗粒大小等几乎不受各向异性影响或无明显相关性，但与迂曲度呈明显负相关。②三维非均质性多孔介质模型中，Rφ相同，其连通孔隙度、孔径受非均质性的影响较小，其值波动幅度不大；迂曲度与非均质性虽无明显相关性，但其值波动幅度明显；多孔介质非均质性与固相颗粒大小有关，进而影响比表面积，比表面积与非均质性呈明显负相关。

表4 孔隙结构参数与各向异性、非均质性的相关性

多孔介质的孔隙尺度各向异性显著影响迂曲度，孔隙分布非均质性显著影响比表面积，进而显著影响多孔介质中流体的流动动态：①主流方向各向异性因子越大，多孔介质中孔隙长轴沿主流方向的取向性更强（孔隙长轴更倾向于朝向主流方向）[21]，流动路径更趋平直，迂曲度显著减小，流体流过多孔介质所消耗的能量减少，多孔介质的渗流能力增强，迂曲度与各向异性的强相关是各向异性影响渗透率的根本原因。②孔隙分布非均质性增强，多孔介质中固相颗粒尺寸越大，孔隙聚集现象更加明显[26]，孔隙间发生相互重叠的概率增加，比表面积将显著减小；通常多孔介质孔隙度相近，而固相颗粒越大或比表面积越小时，流体流过多孔介质所受到的壁面摩擦阻力越小[38]，多孔介质的渗透率将越大。③虽图14显示三维非均质性多孔介质模型的渗透率与比表面积无明显相关性，但总的趋势表现为渗透率随非均质性的增加而增加；这是因为迂曲度也是影响多孔介质渗透率的重要因素[39]，迂曲度与多孔介质中固相颗粒的排列方向和位置（尤其是其中较大的固相颗粒）有关[2]，三维非均质性多孔介质模型中，随着非均质性的增强，生成的固相颗粒间差异越大，分布位置的随机性可能导致迂曲度随机变化，可增大也可减小，导致迂曲度与非均质性无明显相关性；孔隙分布非均质性对多孔介质渗透率的影响是比表面积和迂曲度共同作用的结果。④采用迂曲度与比表面积的乘积分析其与非均质性的关系（见图19）显示，其与非均质性呈明显负相关，说明迂曲度与比表面积共同作用，是非均质性影响多孔介质渗透率的内在原因。

图19 非均质性模型迂曲度、比表面积乘积与非均质性因子的关系

对模拟结果进行数值拟合，得到三维各向异性、非均质性多孔介质模型（未细分Rφ）的渗透率与迂曲度满足幂函数关系（见表5）。

系数a和b与多孔介质孔隙结构有关，且b值与多孔介质的孔径呈负相关，与比表面积呈正相关。拟合结果显示相关系数平方值均接近70%，最高超过80%，精度较高，可用于岩心渗透率的近似估算。

表5 三维多孔介质的渗透率与迂曲度关系式

综上所述，为更好地描述各向异性、非均质性对流体流动的影响，表征非常规储集层的渗透率模型必需包含孔隙尺度各向异性、孔隙分布非均质性这2项反映其微观孔隙结构性质的参数[15]。

4 结论

致密砂岩孔隙形态复杂，孔隙尺度各向异性、孔隙分布非均质性显著，其中主流方向各向异性因子最大，其他方向相近。

各向异性影响多孔介质中孔隙长轴的取向性及流体流动路径，沿各向异性因子大的方向迂曲度小、流体流动消耗能量小，迂曲度与各向异性的强相关是各向异性影响渗透率的根本原因。

非均质性对渗透率的影响表现为迂曲度与比表面积的共同作用，比表面积、迂曲度的乘积与非均质性呈明显负相关，孔隙分布非均质性越强，乘积值越小，渗透率越大。

复杂多孔介质的渗透率与迂曲度满足乘幂关系式，拟合精度较高，可用于岩心渗透率的近似估算。

符号注释：

a——拟合系数，10-3μm2；A——各向异性因子，无因次；b——拟合系数，无因次；cs——声速，m/s；Cs——固相生长核；d——固相颗粒大小，m；D——固相生长速度，f；Dx、Dy、Dz——x，y，z方向的固相生长速度，f；ei——格子速度，m/s；fi（x，t），fj——离散分布函数，kg/m3；feq,j——平衡态分布函数，kg/m3；g——孔隙或子单元块的序号；H——孔隙分布非均质性因子，无因次；hg——第g个孔隙的高度，m；i，j——离散速度方向；k——QSGS算法中固相生长方向；K0——渗透率，10-3μm2；M——转换矩阵；M-1——M的逆矩阵；n——多孔介质划分的子单元块数量，块；nx，ny——一条沿x或y方向穿过多孔介质的直线，其单位长度所遇到的平均孔隙数，个；N——多孔介质中孔隙总数量，个；o，p，q——多孔介质中x，y，z方向的坐标；Pc——粗结构固相生长核概率，f；Pcd——固相生长核概率，f；Pk——不同方向的生长速度，f；Pr——细结构固相生长核概率，f；r——孔隙半径包含格子数，个；Re——雷诺数，无因次；se，sm，sq，sv，sε，sπ——能量矩、三阶矩、能量通量矩、应力张量矩、能量的平方矩、四阶矩的弛豫速率，无因次；S——比表面积，cm2/cm3；S——对角矩阵；t——时间，s；u——流体速度，m/s；u（o，p，q）——多孔介质（o，p，q）处流体速度，m/s；u（o，p，q，t）——t时刻多孔介质（o，p，q）处流体速度，m/s；ux（o，p，q）——多孔介质（o，p，q）处x方向的流体速度，m/s；——x方向平均流体速度，m/s；wg——第g个孔隙的宽度，m；wi——权系数，无因次；x，y，z——空间直角坐标系的3个方向；δt——时间步长，s；δE——模拟结果相对偏差，无因次；δL——单位网格长度，m；∇p——压力梯度，Pa/m；μ——动力黏度，Pa·s；ν——运动黏度，m2/s；ρ——流体密度，kg/m3；τ——迂曲度，无因次；τR——无因次弛豫时间；φ——多孔介质总孔隙度，%；φg——第g个子单元块的孔隙度，%；φeff——连通孔隙度，%；φR——细结构孔隙度，%。