岩石孔隙结构表征与流体输运可视化研究进展

岩石作为一种非均质性的多孔介质材料，其内部复杂无序的微观孔隙结构形态与空间分布对岩石的宏观力学性能(强度、变形特性)、输运特性(流体流动、传热传质)和波动特性等具有极强的控制作用[1]。大量的理论分析、室内实验及现场工程实践表明，诸如石油工程开发，大型水利水电设施建设，地下水开采，隧道和地下硐室开挖，以及地下核废料处理与CO2地质封存等重大工程实践都与岩石特性密切相关。因此，开展多孔岩石材料的微观结构表征有助于深入理解岩石宏观物理现象及作用机制，同时，对于保障工程建设的安全开展具有重要的理论和实践指导意义[2]。

多孔岩石的微观孔隙结构是岩石中流体(包含气体和液体)流动和赋存的场所，微观孔隙结构特征不仅对岩石内部流体的输运特性有重要影响，同时也直接影响介质的宏观物理特性。油气田开发工程中，对流体在岩石中的渗流机理研究是准确评价油气产能和高效开发的理论基础。长期以来，对油气储层流体的渗流规律研究主要通过室内实验和理论分析相结合的方法，利用天然取芯或人造岩芯开展室内流动实验[3]，从而获得可以指导油藏工程开发的各项物理参数。然而，随着低渗致密储层开发程度的加大，此类储集岩复杂的微观孔隙结构、极强的非均质性和多尺度特征使得其许多宏观特性(如毛管压力、孔隙度和渗透率等)同常规的中高渗储层存在显著差异，特别是此类储集岩中的流体流动表现出极强的非线性[4]。因此，加深此类储集岩的微观结构特征和微观孔隙结构对渗流规律的影响机制研究，对此类储层的有效开发具有实际指导意义。

但是，大量的实验结果表明，由于此类储层大量发育纳米级的孔隙，同时受限于室内实验测试技术与测试设备，目前，常规的室内实验测试手段还难以精确表征低渗致密储层的微纳米尺度特征，渗流机理的研究也停留在宏观的描述，且相关的测试技术复杂、实验周期长以及测试成本昂贵等因素极大地制约了低渗致密储层的渗流机理研究[5]。因此，开展基于高分辨率岩芯CT图像的数字模型重建与数值模拟方法研究，有助于微观孔隙结构的精细化定量表征与微观尺度的渗流机理和流体输运特性的分析与评价。通过微观孔隙结构的表征实现渗流通道的精细刻画、研究微观尺度上的渗流机理，为此类低渗致密储层的产能评价和后续开发方式的设计奠定坚实的理论基础。

本文从岩石微观孔隙结构的室内测试与表征、微观可视化渗流物理实验以及基于岩芯CT图像的孔隙尺度模型重建和数值模拟方法(数字岩芯技术)等3方面出发，综述了储层多孔岩石微观孔隙结构实验测试与定量表征技术，多孔岩石输运特性的物理实验与数值模拟方法的研究现状、取得的进展，以期为多孔岩石微尺度孔隙结构表征和输运特性可视化研究提供指导。

岩石孔隙结构是对岩石具有的大量复杂且无序的孔隙和喉道的几何形状、几何尺寸、空间分布及相互连通关系的总体描述。孔隙反映了岩石的储集能力，喉道的大小、形状及连通性则决定了岩石的渗透性[6]。大量的勘探开发实践表明，储集岩的微观孔隙结构直接影响着储层油气的储集与渗流能力，从而控制着油气藏产能的差异性分布[7]。

目前，随着实验测试设备和分析手段的发展，孔隙结构获取与表征的室内测试技术已经日趋完善与丰富，形成了包含孔隙结构特征定性描述、半定量分析与定量表征为一体的分析测试体系。从最初的常规压汞和铸体薄片分析，逐步发展到结合各种先进设备的恒速压汞、核磁共振(Nuclear Magnetic Resonance：NMR)、扫描电镜(Scanning Electronic Microscope：SEM)、背散射电子成像(Back Scattered Electron：BSE)、原子力扫描电镜(Atomic Force-SEM)、聚焦离子束扫描电镜(FIB-SEM)、微—纳米CT扫描成像(Micro-Nano Meter Computed Tomography：μ-CT)，以及利用测井资料进行解释等研究方法(表 1)。理论分析方法也呈现出地学、化学、数学和物理等多学科交叉的特点。

早期研究中，通过常规压汞实验获取进汞压力与进汞体积之间的关系得到孔喉之间的综合毛管力曲线，从而获得孔喉的定量分布，但无法将孔隙和喉道区分开[8]。恒速压汞通过控制进汞速率，以近乎准静态的进汞过程分别表征孔隙和喉道，单独给出孔隙和喉道的毛管力曲线，实现孔隙和喉道尺寸分布的定量表征[9]。但是，有学者研究认为[10]，压汞法这类流体注入的测试方法可能对部分纳米级孔隙的孔隙度测试造成缺失(压力限制)，其次，由于遮挡效应(表面狭窄孔喉可能导致内部较大孔喉流体无法进入)和边界效应(压汞在很低压力下进入狭窄孔喉内)对测试结果的精度会造成一定影响。因此，邹才能等[11]提出结合N2和CO2吸附，实现对部分缺失孔隙参数的测定(如大小和孔径分布)，从而形成了以高压压汞技术结合N2和CO2吸附的新流体注入法，弥补了压汞法存在的不足。核磁共振分析通过测定实验样品中流体的流动来反演介质的微观孔隙结构特征，可用于评价结构特征及流体参数对流动特性的影响[12]。此外，通过测井资料解释也能够有效、快速地分析岩石的孔隙结构特征，推算孔隙结构参数的定量表达式[13-14]。

岩芯铸体薄片是指在真空状态下，将染色的树脂或胶状物加压填充到岩石孔隙中，通过固化切片处理得到能够用于在显微镜下观测的岩石二维薄片，从而实现孔隙特征识别、孔喉连通性分析及矿物组分鉴定[15]。扫描电镜通过电子束对固态样品表面的二次激发获取能够反映标本表面结构的扫描图像(图 1)，可用于岩石表层形貌、矿物组分等相关研究[16]。原子力扫描电镜可以对样品进行原子级分辨率形貌探测，提高了样品形貌和微观结构的成像精度，目前在生物和材料科学领域应用较为广泛[17]。聚焦离子束扫描电镜可以通过离子束对样品进行切片再逐层扫描的操作实现三维结构解析，从而在扫描电镜的基础上极大地提高微观结构的探测精度[18]。

以上提及的各项测试都需要对测试样品进行前处理，完成测试试样的制备，因此，对于样品的原始微观孔隙结构不可避免地会产生一定程度的破坏，得到的测试结果也存在一定失真。微CT扫描成像通过X射线入射(图 2)，获取样品的投影数据，结合图像重构算法完成内部结构的重建，因此是一种真正的无损三维成像测试技术，可广泛应用于岩石介质内部孔喉结构形态、尺寸分布、流体流动和驱替等实验研究中[19-20]。

以上提及的各项测试技术手段根据其测试原理与方法，获取的数据与表征形式有所不同，但是都以获取孔隙结构特征为目的。其中以压汞法和气体吸附法为主的流体注入法主要以获取孔喉特征及孔喉半径等统计特性数据为目的实现孔隙结构的评价与定量表征；通过核磁共振信号结合测井资料构建核磁毛管圧力曲线的反分析法实现孔隙结构的评价与定量表征。这两类方法主要获取孔隙结构的统计性特征分布实现储层孔隙结构的评价与表征，并不能获得直观的孔隙结构图像与岩芯矿物组分分布等特征，难以开展后续的相关测试与研究工作。铸体薄片、扫描电镜和CT扫描等光学辐射法则能获得孔隙结构及矿物组分分布的图像，除了能更加直观地对孔隙的分布进行观测与描述，结合数字图像分析技术还能够实现孔隙结构特征的定量表征。随着实验设备和测试技术理论的发展，除了能够更加定量与精细化地评价和表征岩石微观孔隙结构和矿物组分分布的微结构特征，同时也能获取更高分辨率的微尺度图像，为实现高精度数字模型重建与数值模拟分析奠定基础。

1856年提出的达西定律反映了流体在多孔介质中流动的线性规律[21]，也奠定了渗流力学模型实验的基础。渗流物理模拟实验对油气田开发工程中储层流动的机理分析及开采机理的研究具有重要的理论及实践意义。根据物理模型尺度及研究意义的不同，渗流物理模拟实验主要分为微观尺度的微观可视化渗流模型实验和宏观尺度的室内岩芯及大尺度三维物理模拟试验[22-23]。前者主要用于孔隙空间内流体流动及驱替现象的微观机理分析，后者多用于油藏渗流规律和指导实际工程开发的开采机理研究。

室内岩芯尺度流动分析及驱替实验是目前最为成熟、应用最为广泛的渗流物理模拟技术。室内岩芯流动测试是岩芯分析的一项基础工作，可以对岩芯微观孔隙结构和流动特征进行精确表征与分析[24]。室内岩芯测试的模型以天然柱塞岩芯为主，对于物性较差、取芯成本较高的储层，也可以采用人造岩芯代替[25]。岩芯渗透率是表征储层流动特性的一项基础参数，目前，室内岩芯渗透率测试主要分为稳态法和非稳态法两大类[23]。室内岩芯尺度的流动分析和驱替实验已经广泛地用于岩芯渗透率、相对渗透率和毛管压力等油藏开发基础参数的测定，同时也用于聚合物驱、三元复合驱、泡沫驱和表面活性剂驱等提高采收率技术的室内验证以指导油藏开发[23]。

大尺度的宏观物理模型实验主要是基于相似理论，利用大尺寸的三维物理模型对油藏进行缩放，从而进行井网优化、水驱、化学驱和热采等油藏开发模拟，更加侧重实际油藏条件下的流体饱和度、压力场和温度场的动态变化[26]。大尺度物理模拟实验主要通过大理石组合、透明砂箱堆积及有机玻璃板刻蚀等方法来构建宏观的物理模型，根据模拟对象的不同可以分为流动机理模型和开采机理模型。前者基于理想的规则结构来模拟流体流动规律，后者基于复杂的裂缝网络来模拟井网间的流动特征。大尺度物理模型实验主要用于碳酸盐岩这一类缝洞型油藏的渗流机理研究，这类物理模拟主要通过砂箱堆积、大理石板组合及有机玻璃刻蚀等手段来制备。按模型缝洞的组合形式可以分为流动机理模型和开采机理模型，前者通过规则的裂缝和溶洞模拟油水流动规律，后者通过非均质裂缝网络实现开采机理的研究[27]。由于传统模型制备方法对缝洞形态与尺寸难以实现精确控制，如裂缝开度与倾角、裂隙密度、缝网布置形式、溶洞洞径和溶洞密度等，因此，现有大尺度流动模拟实验大多停留在定性研究的阶段。

但是，不论是岩芯尺度的渗流模型实验，还是大尺寸三维物理模拟系统，由于模型介质的不透明，同时受限于测试技术与现有测试设备，目前还难以在实验中直观地获得模型内部流体流动、流体分布、液—液或固—液界面间相互作用、岩芯结构变化等现象的图像化观测结果。而以上提及的种种现象在岩芯流动或驱替过程中广泛存在，且极大地影响着岩芯流动行为及驱替效率，实时捕捉流体分布与岩芯结构的动态变化，研究其中的微观作用机制，对油藏开发的动态设计具有重要的理论意义。因此，基于真实孔隙结构制备的微观渗流物理模型为研究流体在孔隙空间内的复杂流动行为及渗流过程的实时可视化观测提供了可能。

微观渗流物理模型实验是流体渗流机理分析及驱替现象可视化观测的孔隙尺度实验方法，从微观尺度上揭示了储集岩中的流体渗流规律，为提高采收率研究提供了理论依据和技术指导。典型的微观渗流实验物理模型的制作是利用具有较好连通性的孔隙结构二维图像为模板，利用显影技术将孔隙结构曝光在玻璃薄片上；然后，通过氢氟酸等化学试剂定向腐蚀或激光刻蚀的手段，将孔隙结构的轮廓复刻到玻璃薄片；最后，通过模型盖板的粘接、液体注入与流出通道的制作完成模型的制备，得到具有真实孔隙结构的微观渗流物理模型[22]。另一种方法是通过翻模硅胶等具有较高可塑性的材料，仍然以二维孔隙结构图像为模板，通过制模、翻模、粘结成型等技术制备物理模型[28]。

微观渗流实验的开展需要利用显微镜，通过精密电子仪器完成流动参数(包括流量、压力等)的测定以及镜下流动图像的捕捉(图 3)。微观可视化渗流模型实验的发展最早可以追溯到20世纪50年代[22]。微观渗流模型的发展经历了最早的夹珠模型[29]，过渡到毛管网络模型[30]，然后是与真实孔隙结构特征相似的孔隙网络模型[31]和真实砂岩模型[32]。微观渗流模型已经被证实是一种研究孔隙尺度下单相流体流动机理、多相流体相间作用、固—液界面物化反应和提高采收率机理的可视化研究工具。因此，微观结构特征与自然多孔介质的相似性和良好的光学透明性是微观渗流模型实验技术的重要实施条件[23]。

随着微观渗流物理模型实验技术的发展，众多学者基于此项技术开展了大量微观渗流及微观驱替(如水驱、气驱、聚合物驱和泡沫驱等)提高采收率方法的室内实验(图 4)，系统总结与完善了岩芯尺度下不能具体和精确描述微观渗流规律的问题。同时，基于微观驱替实验中卡断、绕流、指进现象[33-34]的可视化观测，加深了对于剩余油形成与分布、最终采收效果的微观机理认识[35]。如微观水驱油实验中，驱替相液体从大孔周围连通性更好的小孔里流走，从而呈现出“小孔包围大孔”的特殊现象，这也解释了为什么在均质储层开发中出现大面积剩余油的原因；注聚合物提高采收率的微观实验结果表明，聚合物不但能提高波及系数，还能提高洗油系数，这一发现为大规模推广聚合物驱提高采收率的开采技术提供了新的证据。此外，该技术也为储层岩石润湿性、孔隙结构特征和流体性质对采收率的影响机制提供了科学解释，为合理选择开采工艺、实现油气田的有效开发提供了科学依据。

微观渗流物理模型实验能够在形态、尺寸及微观结构方面较好地模拟真实岩芯，实验结果的清晰度也较高，能直观且较为精确地反映孔隙空间中的渗流过程，可以有效地帮助我们理解流体在微观尺度下的流动规律。但是，微观薄片模型的制作与渗流实验操作是开展此项技术的一个难点，更重要的是，二维的平面径向渗流规律和真实的三维空间渗流规律存在本质的不同，孔隙结构的空间连通性使得渗流过程变得更加复杂。因此，如何实现孔隙尺度下流体在三维空间中的可视化观测，对还原流体真实赋存状态与空间流动行为的模拟具有重要意义。鉴于此，研究人员将核磁及超声波探测等技术引入流体流动实验中，以期实现流动过程的准实时监测。

目前，结合低场核磁共振分析或超声波探测技术可以实现流动实验中流场的可视化监测、岩芯内部流体分布及饱和度计算等[36-37]。但这类方法并不能满足岩芯驱替实验中流体分布及岩芯结构变化的实时监测，并且对于弛豫时间接近的两种流体，在将核磁信号转换为图像数据时，由于不同相流体在图像中的对比度低，因此，也难以实现这类混相流体的识别与分割。鉴于此，通过引入CT扫描技术，可以在不改变岩芯外部形态和原始内部结构的条件下，通过对实验装置内不同岩芯部位以及不同驱替阶段的成像，获取驱替实验过程中岩芯内部微观结构、流体流动和分布的动态演化特征[38-39]。结合实验中的压力、流量等监测信息及CT切片扫描技术，可以获取岩芯内部流体饱和度及非均质性的沿程分布特征，从而实现微观驱油效率、剩余油分布特征的定性描述与定量表征[40-43]。

准实时CT扫描实验中，CT扫描也分为几个主要模块，分别为扫描干岩样获取岩芯孔隙结构，扫描饱和地层水的岩芯获取流体初始分布状况；模拟排驱地层水时，同时检测压力降，按一定时间间隔进行CT扫描，直到孔隙流体形态和压力不随时间变化即认为达到油饱和，利用CT图像分析建立初始束缚水饱和度；最后是模拟驱替过程的CT扫描，当岩芯饱和度分布和流体赋存形态不随时间变化，即认为驱替结束，通过CT成像获取剩余油饱和度和流体赋存形态及其分布(图 5，图 6)。

得益于模型可重复使用、表面物化特性可控等优点，微观渗流模型实验可模拟不同的地下环境，再现真实的渗流过程，在认识储层渗流规律和提高采收率研究方面具有其独特的作用。基于准实时CT扫描的渗流过程动态监测将二维的可视化观测扩展到了三维空间，更加符合真实的流动状态，虽然存在模型尺寸与成像精度的矛盾，但也已经迈出了渗流过程实时可视化监测的一大步。随着实验测试硬件和技术的发展，微观渗流模型实验的适用条件及研究范围也将日益广泛。

然而，现阶段的微观渗流物理模拟技术受渗流实验装置、成像设备和扫描技术的限制，目前可查阅的、结合CT技术的微观渗流实验研究，准确地说是一种阶段性、准实时的可视化实验技术。因此，受限于目前的技术水平，微观可视化渗流实验中还难以实现真正的全场实时监测与成像，此外，在这类技术中施加高温和高压环境，模拟真实岩石的地下赋存状态也是一个很大的挑战。因此，基于微CT图像的孔隙尺度建模和计算机模拟的数字岩芯技术就成为了解决室内实验中硬件限制的一个突破口，并且这项技术的优势也越来越明显，也是近年来地球物理学研究的一个热点领域。

近年来，基于数字岩芯技术的岩石微观孔隙结构表征与输运特性研究发展尤为迅速[45]。数字岩芯技术基于虚拟的数字化模型和计算平台，通过微CT成像、建模和计算机模拟揭示传统实验无法展现的岩石内部物理现象(图 7)，从微观角度解释宏观现象，实现油藏的精确描述。数字岩芯技术易于实现多物理场及多尺度效应的耦合，同样也可以实现模型参数及实验微观影响因素的定量控制。相较于传统的室内实验方法，数字岩芯技术的发展不仅丰富了传统室内实验测试的研究手段，因其具有实验条件可控、实验效率高、模型可重复利用且实验成本较低等优点也成为了地球物理领域的研究热点[46-48]。数字岩芯技术已经在复杂多孔介质输运特性预测的研究中取得了卓有意义的成果，同样在岩石物理参数，如声学、电学和热传导等特性的研究中具有良好的应用潜力[49-50]。

作为数字岩芯技术的基础工作，构建能够反映真实岩芯微尺度结构特征的数学计算模型是保证后续模拟成功开展及结果真实可靠的保障。根据模型构建方法的不同，微尺度模型可以分为非重建模型和重建模型两大类。前者基于规则的几何单元来表征复杂多孔介质中的孔隙结构，这类模型也称为理想化模型；后者基于真实岩芯的孔隙结构图像数据，根据不同的重建算法完成模型构建，这类模型也称为真实模型[51]。

非重建模型主要通过给定规则的几何结构及尺寸参数来构建多孔介质模型。最早的非重建模型可以追溯到20世纪50年代，Fatt等通过在二维规则晶格中随机给定格子半径，利用这种简单的网络模型来研究地下水的毛管压力和渗透率[52-53]。受限于计算能力，最初的晶格模型利用电阻器物理模型进行模拟，流场等效于电场，水力传导系数等效于电传导率。但有学者指出，这类模型无法表征多孔介质的空间连通性[54]。

平行毛细管束模型的概念最早于1949年由Purcell提出[55]，1975年Dullien将该模型应用于单相流体的流动分析[56]。该类模型将岩石孔隙等效于一系列不同直径的平行毛细管束，可以表征孔隙的空间连通性，由于采用等直径的管，因此无法反映孔隙间的交叉连通性。为了进一步改进毛细管束模型，能够更加精确地反映多孔介质拓扑结构特征，通过改变平行毛管的排列分布，引入含有分叉的毛管束，以及与水阻力、毛管力相关的函数来控制模型内圆管分布[57]；或者结合分形理论，构建含有变直径弯曲毛细管的改进毛细管束模型，该类模型多采用考虑毛管压力的泊肃叶定律研究多孔介质中复杂流体性质和流动数学模型的推导[58-59]。

在岩石微观孔隙拓扑结构的研究方面，堆积球模型的提出极大地改善了模型的拓扑结构特征表征，并再现了多孔介质的空间水力连通性[60]。该类模型综合考虑了地质演化过程中岩石的沉降、压实和胶结作用，因为球体的中心位置已知，可以通过Voronoi Tessellation(曲面细分)方法提取堆积球模型中的孔隙网络结构。该类模型可用于多孔介质中单相和多相流体的渗透率、相对渗透率和毛管压力等输运特性的研究[61]。

另一类非重建多孔介质模型利用一系列规则的二维几何体(如圆形、正方形和三角形等)来代替多孔介质的骨架颗粒，从而形成孔隙间相互连通的多孔介质模型[62]。此类模型可以利用成熟的有限元软件，基于求解纳维尔—斯托克斯(N-S)方程来模拟孔隙中的流体运移及骨架对渗透率的影响。将此类模型扩展到三维空间，利用同样规则的几何体来表征孔隙和喉道(球体代表孔隙，圆柱或棱柱代表喉道)，通过改变与孔隙连接的几何体数量，就可以得到不同配位数的理想化孔隙网络模型[63]。此类理想化的孔隙网络模型，虽然较之以往的结构有了很大改进，但仍然难以反映真实多孔介质中孔隙的拓扑结构特征和实现空间位置的精确定位。

上述各种多孔介质理想化模型如图 8所示，从本质上来说，非重建模型都是晶格模型和毛细管束模型的演变。此类理想化模型虽然能够在一定程度上反映岩石介质内部的结构特征，但究其模型的构建原理，并不能反映岩石真实的孔隙结构拓扑特征，孔隙与喉道的连接及空间分布也与真实结构相差较远。因此，这类模型主要用于理想化的流体输运特征分析，推导适用于宏观渗流分析的经验公式。此外，后续基于重建模型的数值分析计算也往往采用这类理想化模型推导的数学表达式。

微尺度重建模型通过获取天然岩芯的微观孔隙结构图像，通过不同的模型重建算法实现数字图像到可用于数值模拟计算的数学模型构建。根据模型重建原理的不同，目前，微尺度重建模型的构建方法主要有两大类：数值重建法和物理实验法[64-65]。前者一般基于二维的岩芯切片图像分析，利用数学方法构建具有统计学相似性的三维孔隙尺度重建模型，数值重建法一般需要很少的重建数据，具有速度快、成本低等优点；后者基于CT扫描三维成像技术，借助计算机图像处理构建可以直接反映岩芯微观孔隙结构特征的三维重建模型，相较于数值重建方法，此类模型能够反映天然岩芯真实的拓扑结构特征，且能精确地定位孔隙和喉道的空间位置，并表征孔隙和喉道的分布。

目前，数值重建法主要分为过程法和统计法两大类，如图 9所示，前者主要用于模拟成岩过程简单的岩石类型，后者通过统计图像的孔隙度信息及参数自相关函数，通过反复迭代的方式来实现三维模型的重建。常用的数值重建方法及其特点如表 2所示。尽管基于二维图像的数值重建方法需要的数据较少，且模型构建速度也较快，但此类重建模型更多地是表征岩芯微观结构随机统计特征，很难反映孔隙结构真实的空间分布特征。

不同于数值重建的方法，基于三维CT图像提取的孔隙网络重建模型可以真实地再现岩石微观孔隙结构的拓扑特征及孔隙、喉道的空间分布。此类重建模型根据构建算法与模型提取方式的不同又可分为等效孔隙网络模型和孔隙网格模型[51]。

等效孔隙网络模型基于真实的三维微观孔隙结构图像提取孔隙网络结构，并通过理想的规则几何体单元来代替孔隙和喉道。根据提取算法的不同，又可分为中轴线算法、最大球算法和遗传算法[67-69]等。上述提取算法中又以最大球算法应用最为广泛，Blunt等基于最大球算法开发的等效孔隙网络模型提取代码及相应的两相流分析代码被众多学者采用，并在绝对渗透率和相对渗透率预测[70]、模型润湿性与残余油分布研究等方面进行了大量工作[71]。

上述模型也存在不足之处，如在原始图像二值化处理过程中，部分中等RGB值的区域难以界定为骨架还是孔隙，这部分区域往往是由于分辨率的原因而没有识别出的微孔隙或高密度矿物部分，因此，单纯地认为是骨架或者孔隙都不准确。因此，有学者提出了双重孔隙网络模型[72]，通过图像灰度值实现骨架、微孔和大孔的三相分割，提高了模型的预测精度。等效孔隙网络模型从孔喉数量、尺寸和分布上再现了真实多孔介质的微观结构，但采用理想的规则几何体来进行孔喉的表征，因此，在孔隙拓扑结构的表现上也存在不足，虽然引入了形状因子的概念来修正理想假设带来的误差[73]，但仍然无法完全还原孔隙结构的真实几何形状和空间分布。基于等效孔隙网络模型开展的研究取得了很多成果，但这类模型由于自身的局限仅仅涉及了流动分析，而对于宏观渗流过程中存在的固体变形和传热传质等问题考虑不足。

孔隙网格模型利用图像处理和网格划分技术进行微观孔隙模型的重建，基于有限元和有限体积等数值计算方法开展数值模拟，随着计算机硬件水平的提高，这类模型在微尺度建模和数值分析领域具有广阔的应用前景。二维有限元网格模型可以通过提取孔隙轮廓构建几何模型，再利用商用软件实现网格剖分，该方法可以再现孔隙的真实几何结构[74]。这类模型在从二维向三维扩展时，由于数据处理量的激增，将会占用大量的计算机内存(如400体素的重构模型有1~GB左右)。受限于计算机的处理能力，在PC机上采用一般的有限元软件对1~mm左右物理尺寸的模型进行网格划分都有很大难度，增大模型尺寸则几乎不可能完成[75]。因此，为了解决这一个问题，有学者将商业CT图像处理软件(如Amira/Avizo，Mimics和Simpleware等在医学和材料领域广泛使用的图像处理与模型重建软件)引入到多孔介质微观模型的重建中，该方法直接跳过几何模型的构建，在提取的孔隙网络实体表面生成面网格，在此基础上，完成体积单元网格的构建[76]，此类方法生成的网格为三角形或四面体，称为非结构化网格模型。这类方法在网格划分时，为了保证体积单元网格生成的成功率，在微小孔隙和尖端区域进行了优化处理，采用局部网格收缩，膨胀或剔除质量较差的网格等方式，改变了孔隙特征，但生成的网格质量难以保证，因此多用于网格质量要求不高的固体变形数值仿真。

针对以上模型重建方法存在的不足，提出一种结构化网格模型的构建方法：基于网格单元体与图像像素体在空间位置的一一对应关系，以具有相同像素尺寸的网格单元体替换图像像素，构建四边形或六面体的结构化网格模型。结构化网格模型提供了一种稳定、高效的基于微CT图像的模型重建方法，解决了传统重建模型网格质量难以控制、模型拓扑结构失真等问题，并且可再现天然岩石复杂的微结构特征，具有良好的网格质量，可同时满足渗流和力学变形的数值模拟分析，模型构建的具体细节可参考笔者以往的研究工作[77]。

以上提及的关于多孔介质重建模型的物理实验法分类如图 10所示。基于岩石微CT图像的模型重建为岩石微观物性参数预测、宏观变形及流体输运特性的模拟提供了一种快速、稳定的解决方案。相较于传统的室内测试手段，微尺度重建模型具有实验效率高、成本花费低、模型可重复及模拟中部分实验及环境参数可控等诸多优点。同样，如何构建可以真实反映岩石微观孔隙结构特征的重建模型也是这类技术实施的关键，如何提高基于微尺度重建模型的数值模拟计算精度和收敛速率也关乎计算结果的可靠度。因此，在后续研究中，重建模型除了要能够真实反映岩石微观拓扑结构特征、优化模型特征参数(如网格单元数)，提升数值计算效率也是这类技术的一个重点关注方向。

多孔介质中流体输运性质的预测是进行流体宏观流动特性分析的基础。根据不同的表征思想，多孔介质内流体输运方程可以分为分子动力学(Molecular Dynamics，MD)方法[78]、格子玻尔兹曼方法(Lattice Boltzmann Method，LBM)[79]及求解各种偏微分方程的计算流体动力学(Computational Fluid Dynamics，CFD)方法等[80]。

MD方法通过统计模型中所有分子的运动规律来表征流体整体的输运性能，该方法可以精确描述模型中流体粒子的分布情况，进而得到模型中流体的吸附、扩散等性质。近年来，该方法被广泛应用于吸附解析、扩散、微流体和微传热等领域，取得了一系列研究成果[81]。特别是随着页岩油气和致密油等非常规油气的勘探开发，这类储层中的最小孔隙已经达到了纳米尺度，常规的连续性假设在分析这一类储层中的流体流动机理时已经不再适用，转而采用平均分子自由程和孔隙半径比值的克努森数(Knudsen)来判定纳米孔隙中分子的输运模式。纳米孔隙中流体的运移模式主要有分子与分子间碰撞产生的黏性流、分子与壁面间碰撞产生的克努森扩散及壁面吸附分子沿壁面蠕动产生的表面扩散[82]。图 11所示为纳米孔道内流体运移的分子动力学模式示意图。

分子动力学模型与真实多孔介质的无序孔隙结构差距较大，采用的数学方程也相对简单，其核心思想是基于经典的牛顿力学定律，通过对每一个粒子运动方程的求解来获得粒子位置的变化，基本方程为

式中：mi—粒子i的质量，g/mol；

ri—粒子i的位置矢量；

V(rij)—粒子i和j的相互作用势能函数；

rij—离子i和j之间的相对位置；

t—时间，ps(皮秒，1 ps=10-9 s)。

目前，MD方法主要用于纳米尺度下流体输运特性的机理分析，也以单管构建的简单模型为主，在复杂多孔介质模型中的MD方法应用还少有研究报道。

LBM方法是一种将宏观流体力学与微观分子动力学联系起来的介观理论，通过求解离散介观尺度Boltzmann方程来再现微纳尺度的流体流动问题，在低马赫数条件下可以转化为Stokes方程[83]。由于LBM方法的简便性与高计算效率，常用于研究多孔介质中孔隙尺度流体流动问题。例如，该方法可以直接在分割后的CT图像的体素格点上执行，并且不需要生成网格，因此非常适合于具有复杂几何边界的多孔介质中的流体流动模拟[84]。根据弛豫时间的不同，LBM模型的表现形式也有所不同，最常用的是单弛豫时间模型(LBGK)[85, b86]和多弛豫时间模型(MRT)[87]。

LBGK模型是目前应用最广泛的LBM模型，不但极大地减少了模型的计算工作量，并且通过一定的转化还可以导出正确的Navier—Stokes(N—S)方程。其中又以Qian等提出的DnQb模型最具有代表性、应用最为广泛[85](n代表空间维数，b代表离散速度数目)，其离散速度演化方程为

平衡态分布函数为

式中：ci—离散格子速度；

τ—松弛时间因子；

u—速度矢量；

u—宏观速度，m/s；

ωi—权系数(与离散速度方向的矢量长度有关)；

ρ—密度，g/cm3；

${c_{\text{s}}} = \sqrt {RT} $，与声速有关；

R，T—与空间维度相关的常数。

MRT模型与LBGK模型的最主要区别在于碰撞过程使用多个松弛时间，可以表示为

式中：$\mathbf{\Lambda}^{-}$—碰撞矩阵；

i = 1, 2, ..., b。

该方程描述了离散分布函数

在空间V = Rb中的时空演化过程，空间V成为分布函数空间或者离散速度空间。关于多弛豫时间模型MRT的基本方程及相关程序算例可以参考文献[88]。

LBGK模型中常用的DnQb模型主要有一维的D1Q3和D1Q5、二维的D2Q7和D2Q9及三维的D3Q15和D3Q19模型。图 12所示分别为二维的D2Q9和三维的D3Q19离散速度网格结构图。

利用LBM计算数字岩芯的渗透率时，离散速度演化方程(式(2))、平衡态分布函数(式(3))、宏观密度$\rho = \sum_{i}f_{i} \left (x, t\right )$与宏观速度$\rho u = \sum_{i}f_{i} \left (x, t\right )e_{i}$等共同构成了多孔介质中流体流动的LBM迭代模型。实际计算中，通过选定流体渗流方向，出口和入口施加一定压差驱动，其余四面以一层骨架格点封装，孔隙与骨架之间往往采用无滑移反弹边界条件，从而完成计算前处理和模型边界条件的施加。

CFD方法基于连续介质假设，可以将流体看作充满整个流场的连续介质，从而可以定义流场空间每一点上的流体密度、速度、温度和压力等。连续介质假设是流体力学的一个基本假设，可以看作一种对流体的近似处理，流体的运动满足质量守恒、动量守恒和能量守恒定律，在数学上可以用一组偏微分方程描述[88]

式中：

e—单位质量的内能，J；

σ—应力张量，Pa；

q—由热传导和热辐射引起的热通量，J/s。

由于式(6)并不封闭，因此需要引入应力张量和变形率张量、热通量与温度梯度间的关系，以及关联热力学性质的状态方程，并且在Stokes假设之下，该方程才完整。基于以上假设，该偏微分方程组又称为N—S方程。如果忽略流体的黏性项，即对于理想流体而言，该方程可以简化为Euler方程。目前，基于连续性假设的CFD方法是流体力学领域发展最为成熟、应用最为广泛的一种方法。目前，常用的数值求解软件，如Fluent、Comsol、ANSYS、Abaqus、CFX等，均是基于N—S方程组的求解。以Euler和N—S方程组为控制方程，通过各种离散方法进行求解。常见的模拟方法包括有限差分(Finite Difference Method，FDM)[89]、有限容积(Finite Volume Method，FVM)[90]、有限单元(Finite Element Method，FEM)[74]、有限分析(Finite Analytical Method，FAM)[91]和边界元法(Boundary Element Method，BEM)[92]等。

三种不同尺度的流体控制方程可以看作是流体同一物理现象的不同表征形式，三者之间可以通过一定的方法进行等价转换。通过对分子动力学模型中流体分子的统计平均可以得到Boltzmann方法中的流体控制方程；通过Chapman-Enskog展开及Taylor级数展开，可以将Boltzmann方程回归到Euler方程或Navier-Stokes方程组[51]。

上文提及的各种孔隙网络模型及数值模拟方法就是为了构建一种数字化的计算平台，在微观基础上更快、更直观地研究和展示多孔介质中复杂的流动机理和流动形态。孔隙级流动模拟可以得到诸如毛细管力、相对渗透率等反映流体流动特征的基础参数，通过与室内实验测试资料的对比可以验证模型构建的可靠性和计算结果的准确性。

孔隙网络模型利用规则几何体组成的空间网络结构来表征多孔介质复杂的孔隙空间，通过模拟流体在这类简单几何体中的流动过程可以近似地再现多孔介质中的流动。孔隙网络模型的基本组成单元———孔隙和喉道，由一系列等截面的柱体(三角形，正方形或圆形截面)构成，因此，整个网络模型中的流动模拟可以看作大量相互连接的毛细管内流体流动的分析。为了修正孔隙网络模型与真实多孔介质的误差，引入形状因子的概念对孔喉截面进行描述[73]。形状因子不仅在孔隙网络模型的抽提中有重要作用，并且对于模型流动分析也具有一定影响。此类孔隙网络模型为准静态模型，即毛管力的作用远大于黏滞力的作用，毛细管内单相牛顿流体流动的分析，在不可压缩和流体黏度为常数的假设下，利用连续性方程和N—S方程可以求得椭圆形Possion方程来描述毛细管内的流动

截面为圆形时，可以利用Poiseuille方程实现流动过程的表征和求解

式中：v—流体速度，m/s；

μ—黏度，Pa·s；

f—单位质量体积力，N/cm3，此为常数；

Q—体积流量，cm3/s；

Δp—沿毛细管压力降，Pa；

r—毛管半径，μm；

L—毛细管长度，μm。

对于其他形状，特别是任意三角形的情况，则无法用解析式直接进行表征。

利用式(7)和式(8)可以实现圆形截面毛细管传导率的求解，通过进一步无因次化计算，其无因次传导率与毛细管半径无关，是圆形截面形状因子的1/2。对于截面为正方形和三角形的传导率计算可以参考文献[93]。

在模拟驱替过程的传导率计算时，圆形截面的计算过程相对简单，然而，对于三角形或正方形截面，由于毛管夹角的存在，不同相液体在毛细管中分布形态复杂，使得多相流情况下毛细管传导率的计算变得更加复杂，如图 13所示。对于这类毛细管中间部位的计算仍然采用前述单相流体的方法，而夹角处的计算，目前一般采用数值模拟方法实现[94]。

对于孔隙级网络采用准静态驱替模型实现孔隙尺度驱替的模拟，驱替过程中由于模型润湿性的转变，使得驱替形式变得复杂，依据不同的驱替形式，又可分为活塞式驱替、被驱替相流体卡断及孔隙充填等[95]。在准静态驱替模拟中，驱替相流体能否进入孔隙和喉道，完全取决于毛细管的入口压力(驱替压差大于毛管入口压力，驱替相才能进入孔隙和喉道)。而孔隙(或孔喉)间的流体交换产生的压力降由前文所述的传导率控制，孔隙和孔喉的毛管入口压力则由孔隙中流体分布状态和驱替过程决定[96]。对于圆形截面的孔隙，在活塞式油驱水下，孔隙的毛管入口压力的表达式为

式中：σow—油水间界面张力，mN/m；

θowr—油水接触角的后退角，(°)。

对于多边形截面的孔隙，由于油水界面不仅存在于驱替前缘，同时存在于角落位置，因此，Oren等[60]基于相界面力平衡原则，推导并改进了此类孔隙的毛管入口压力

式中：G—形状因子；

Fd—无因次校正因子；

β—多边形内角半角，(°)；且β1 < β2 < β3。

对于出现被驱替相卡断时的毛管入口压力的计算存在两种不同情况，一是油水两相界面同时向孔隙中央膨胀发生的卡断，二是最小角落的界面向最大角落膨胀发生的卡断，两者的毛管入口压力分别由式(11)和式(12)给出

式中：θa—油水接触角的前进角，(°)；

θh, 3—最大角落处油水接触角，(°)。

发生孔隙充填时，由于难以准确计算发生充填的相邻孔喉的空间位置，对于配位数为N的孔隙体有N-1种充填可能性。因此，入口毛管力的计算通常采用Blunt于1988年提出的参数化计算模型[97]

式中：

n—发生充填的孔喉数；

Ai—无因次参数，与渗透率满足${A_2}-{A_n} = \frac{{0.03}}{{\sqrt K }}$的关系；

xi—0~1的随机数。

在利用孔隙网络模型模拟流体在多孔介质中的流动过程时，毛管入口压力决定了流体在模型中的流动顺序。毛管入口压力不仅受控于孔喉截面类型，还与孔喉中流体分布形态、驱替过程和驱替形式有关。因此，为了使孔隙网络模型的计算结果(毛细管压力和相对渗透率)更接近实验结果，除了要在孔隙尺度上反映流体的分布形态，还应尽可能地考虑各种微观驱替过程。

综上所述，利用孔隙网络模型进行流动模拟的基本思路可以概括为：初始模型饱和水，模型为亲水性；进行油驱水，模拟初始原油运移，建立束缚水饱和度；油驱水结束，模型壁面润湿性发生不同程度的变换，通过修改油水接触角来模拟这一物理变化过程；然后进行水驱油，模拟注水开发过程。驱替与吸吮过程中均可以计算毛细管压力和相对渗透率等参数，并可以进行多轮次、往复地驱替和吸吮过程模拟，为微观渗流现象研究提供了一种快速、便捷的计算平台。

基于孔隙尺度重建模型的数字岩芯分析技术为储层岩石微观结构的精细化表征和流动特性的可视化研究提供了一种准确、高效的计算手段。并且同一个数字模型可开展不同计算条件下的多次平行计算，为大规模的多场、多相、多尺度和多影响因素下的并行计算与分析提供了可能。该方法基于虚拟的数字计算平台，模拟的实验环境定量可控，并且可以模拟多重物理影响因素的响应机制，深入探讨微观层面的相互作用关系，实现常规物理实验难以捕捉的一些微观现象与实验参数的表征与分析。相较于传统的室内实验分析，数字岩石物理技术在实验测试周期、原始数据采集、数据分析深度与数据利用广度等方面的优势越来越明显，已经逐渐成为地球物理和岩石力学等领域的研究热点。

(1) 储层岩石宏观物理特性受控于微观孔隙结构特征，因此，孔隙结构的精细表征对于岩石物性的预测与流体输运特性的分析具有重要意义。目前，各类测试手段都以获取岩石的微观孔隙结构特征为目的，可以实现微观孔隙结构的定性、半定量或定量描述与表征。但各个单一的测试手段均难以完全反映孔隙结构的特征，因此，在大力发展、升级测试设备的基础上，也应该开展多技术联合测试的方法，以期实现岩石微观孔隙结构的高效、低成本精确表征。

(2) 微观渗流物理模拟实验从孔隙尺度上为分析流体在孔隙空间中的流动行为及微观力学机制提供了可能。随着此类技术的发展，微观渗流模型已经被证实不仅是一种研究孔隙尺度下单相流体流动机理的可视化研究工具，在多相流体相间作用、驱替波及效率和固—液界面物化反应，以及卡断、绕流和剩余油形成及分布等复杂流动现象的形成机理分析中也起到了重要作用。准实时CT扫描成像的引入，为微观渗流物理模拟向三维可视化发展迈出了一步，将为深入认识储层渗流规律和提高采收率研究方面提供更多助力。

(3) 基于孔隙尺度重建模型的数字岩芯技术是近年来发展迅速的一种孔隙结构表征和岩石物理参数模拟预测方法，是室内岩芯分析的高效替代手段之一。作为数字岩芯分析技术的基础，构建能够反映真实岩芯孔喉形态及分布的重建模型尤为重要。微CT扫描成像作为非接触式三维模型重建技术已经广泛应用于岩芯微观孔隙结构的模型重建中，但受成像硬件设备的限制和扫描岩芯尺寸与成像精度之间的矛盾，因此，构建能反映岩石宏观物理特性的孔隙结构重建模型是开展数字岩芯分析的重点，也是数值模拟计算结果可靠性的保证。此外，随着计算能力的提升，开展多尺度数字岩芯构建，进行跨尺度数值模拟分析也是这项技术的发展趋势。