基于孔隙网络模型的气液渗流特性.pdf

基于孔隙网络模型的气液渗流特性美合日阿依穆太力普1，李涛1，顾文波1，蔡武2，薛庆楠1，敬嘉珩1，魏景阳1，王骞1(1.新疆大学电气工程学院,新疆乌鲁木齐830017；2.中国矿业大学煤炭资源与安全开采国家重点实验室,江苏徐州221116)摘要：CO2地质封存是实现大量 CO2减排的有效途径之一，而阐明 CO2在岩石孔隙中的运移规律是 CO2长期安全地质封存的关键。基于岩心微观结构重建的孔隙网络模型，不仅可以反映真实岩心孔喉分布规律，又能体现孔隙空间的分布及其发育特征，在多孔介质多相流研究领域具有十分重要的意义。选用储层较为常见的 Berea 砂岩，基于微焦点 X 射线计算机断层扫描(micro-CT)技术获得二维 CT 图像，在 ImageJ 软件中通过尺寸裁剪、降噪滤波以及阈值分割等一系列图像处理过程，重建了岩心样品数字图像，并在 Matlab 软件中借助提取函数实现了二值化；然后，采用最大球算法，通过搜索最大球、建立最大球连通关系、孔喉的识别和参数计算 3 个主要步骤，提取了实验岩心样品中的孔喉拓扑等效网络，并运用 Amira 软件进行可视化，将孔隙和吼道等效为球杆模型。通过对比不同图像处理与孔隙网络提取方法得到的岩心孔隙度、配位数、孔喉几何尺寸与形状因子等结构参数，确定了能够较好反映真实岩心孔隙空间特征的数字岩心重构及孔隙网络模型提取方法，即模型 Pic7：滤波处理时将 Mean 值设置为 1.5，Median 值设置为 2；阈值分割选择 Otsu 算法且应用最大球法提取数字岩心时将最小孔径设置为 1。确定了 Berea 岩心孔隙结构参数：配位数分布在 30 以内，峰值在 5 左右；孔隙半径分布在 80m 内，集中在 20m 左右；吼道半径分布在 60m 内，集中在 10m 左右；孔隙形状因子和吼道形状因子分布在 0.07 以内，孔隙形状因子峰值在 0.03 左右，吼道形状因子峰值在 0.035 左右。最后，借助帝国理工大学开发的孔隙网络两相渗流模拟程序，用 0.103mol/kgNaCl 溶液表征咸水，在 50、12.4MPa 的储层条件下，开展了基于孔隙网络模型的多孔介质内 CO2驱替咸水的渗流数值模拟，通过分析驱替和吸入过程中毛细管压力曲线和相对渗透率曲线，阐明了孔隙网络提取方法对多孔介质多相流过程的影响，再次验证了孔隙网络模型建模方法的可靠性以及所得建模参数的准确性。关键词：孔隙网络模型；最大球算法；多孔介质；数字岩心；CO2地质封存中图分类号：P618.13文献标志码：A文章编号：02539993(2023)07280211Seepage characteristic of gas water based on pore network modelMEIHERIAYIMutailipu1,LITao1,GUWenbo1,CAIWu2,XUEQingnan1,JINGJiaheng1,WEIJingyang1,WANGQian1(1.School of Electrical Engineering,Xinjiang University,Urumqi830017,China;2.State Key Laboratory of Coal Resources and Safe Mining,China Uni-versity of Mining and Technology,Xuzhou221116,China)Abstract:IthasbeenfoundintherecentstudiesthattheCO2geologicalstorageisoneoftheefficientmethodstoachieve收稿日期：20221020修回日期：20230206责任编辑：王晓珍DOI：10.13225/ki.jccs.CN22.1512基金项目：新疆维吾尔自治区自然科学基金青年基金资助项目(2021D01C089)；新疆维吾尔自治区重点研发资助项目(2022B01033-2)；中央引导地方科技发展专项资助项目(ZYYD2022C16)作者简介：美合日阿依穆太力普(1989)，女，新疆喀什人，维吾尔族，副教授，硕士生导师，博士。E-mail：引用格式：美合日阿依穆太力普，李涛，顾文波，等.基于孔隙网络模型的气液渗流特性J.煤炭学报，2023，48(7)：28022812.MEIHERIAYIMutailipu，LITao，GUWenbo，etal.Seepagecharacteristicofgaswaterbasedonporenet-workmodelJ.JournalofChinaCoalSociety，2023，48(7)：28022812.第48卷第7期煤炭学报Vol.48No.72023年7月JOURNALOFCHINACOALSOCIETYJul.2023alargeamountofCO2reduction,whiletheunderstandingofthetransportlawofCO2intherockporespaceisfundament-altotheflexibilityofasafeandlong-termCO2geologicalstoragescheme.Theporenetworkmodelbasedonthecoremi-cro-structurereconstructioncannotonlyreflecttherealcoreporeandthroatdistributionlaw,butalsoreflectthedistribu-tionofporespaceanditsdevelopmentcharacteristics,whichisofimportanceinthefieldofmulti-phaseflowresearchinporousmedia.TheBereasandstone,whichismorecommoninthereservoir,wasselectedtoobtaintwo-dimensionalCTimagesbasedonthemicro-focusX-raycomputedtomography(micro-CT)technology,andthedigitalimagesofcoresampleswerereconstructedinImageJsoftwarethroughaseriesofimageprocessingprocessessuchassizecropping,noisereductionfilteringandthresholdsegmentation,andbinarizedintheMatlabsoftwareviaextractionfunctions.Then,viathemaximumspherealgorithm,thepore-throattopologicalequivalentnetworkintheexperimentalcoresampleswasextrac-tedthroughthreemainsteps,includingsearchingthemaximumsphere,establishingthemaximumsphereconnectivity,identificationofporesandthroatsandparametercalculation,andvisualizedbyusingAmirasoftwaretoequatetheporesandthroatstoasphericalrodmodels.Bycomparingthestructuralparametersofcoreporosity,coordinationnumber,poreandthroatgeometryandshapefactorobtainedviadifferentimageprocessingandporenetworkextractionmethods,adi-gitalcorereconstructionandporenetworkmodelextractionmethodthatcanbetterreflectthespatialcharacteristicsofrealcoreporeswasdetermined,suchasmodel7:themeanvalueissetto1.5andmedianvalueissetto2duringfiltering.TheOtsualgorithmwasselectedforthresholdsegmentationandtheminimumaperturewassetto1whenthedigitalcorewasextractedbythemaximumspheremethod.TheporestructureparametersofBereacoreweremeasuredasfollows:theco-ordinationnumberwasdistributedwithin30andthepeakvaluewasaround5.Theporeradiuswasdistributedwithin80mandconcentratedaround20m.Thethroatradiuswasdistributedwithin60mandconcentratedaround10m.Theporeshapefactorandthroatshapefactorweredistributedwithin0.07,thepeakofporeshapefactorwasaround0.03,andthepeakofthroatshapefactorwasaround0.035.Finally,basedontheporenetworktwo-phaseseepagesimulationprogramdevelopedbyImperialCollegeofTechnology,theseepagesimulationofCO2-replacedbrineinporousmediawascarriedoutunderthereservoirconditionsof50and12.4MPabycharacterizingbrinewith0.103mol/kgNaClsolution.Then,thecapillarypressurecurveandrelativepermeabilitycurveduringthedrainageandimbibitioncyclewereanalyzedindetail.Theinfluenceoftheporenetworkextractionmethodonthemulti-phaseflowprocessinporousmediawaselu-cidated,andthereliabilityoftheporenetworkmodelingmethodandtheaccuracyoftheobtainedmodelingparametersweredemonstratedagain.Key words:porenetworkmodel；maximumballalgorithm；porousmedia；digitalcore；CO2geologicalstorage在“碳中和”目标下，发展二氧化碳捕集、利用与封存(CCUS)技术不仅是实现化石能源低碳化利用的重要手段，还是保障我国能源安全的战略选择。最常见的 CCUS 类型有 CO2提高原油采收率技术(CO2-EOR)以及深部咸水层 CO2地质封存技术。由于CO2地质利用与封存是多孔介质中伴随相变传热、气液渗流、多工况多因素耦合的复杂过程，导致其在面向商业化推广时仍面临成本高、见效慢等瓶颈。因此，对多孔介质多相流的微观机理研究，包括孔隙空间的发育规模、空间分布等至关重要。近年来，研究者为定量描述微观尺度渗流过程，基于数字岩心，采用格子 Boltzmann 方法或孔隙网络模型进行了孔隙尺度下的多相流动模拟研究1。相比于传统的岩石物理实验，数字岩心技术具有速度快、成本低、精度高、可重复等优点2，能够在微观尺度上定量描述各种因素对宏观物理特性的影响，是一种实际岩心数字化或实际岩心模拟技术。数字岩心建模方法主要有物理实验法和数值重建法 2 类3。数值重建法是基于少量的二维薄片图像，利用其包含的信息，通过相应模拟方法重建三维数字岩心，包括随机模拟法和沉积岩形成过程模拟法4。不同的是，物理实验法需利用先进光学实验仪器对岩心样品拍摄或扫描以获取大量的岩心二维图片，然后用建模软件把二维图片叠加重构成三维数字岩心5-6，由于其周期长且不具重复性，不利于岩心的研究。最常用的图像获取方法有序列切片成像法(SSTM)、激光扫描共聚焦显微镜法(LSCM)和 X 射线 CT 扫描法(micro-CT)6。SSTM 的缺点是在切割和抛光时需要大量的时间以及该过程样品岩心会遭到破坏，而且相邻切片间有一定的厚度间隔，通过一定的算法进行推测无法准确了解切片间孔隙的连通情况7；LSCM 的缺点是不适用于大孔喉的岩样(穿透第7期美合日阿依穆太力普等：基于孔隙网络模型的气液渗流特性2803深度一般不到 1mm)，且孤立孔隙因为无法加入用于增强反射强度的含有荧光剂的环氧树脂而不能被扫描到，导致所建模型准确度较低6。而 micro-CT 具有快速、准确、无损的优点8，故笔者采用 micro-CT 技术提取典型砂岩 Berea 岩心图像。进而建立孔隙网络模型，分析岩心内部孔隙空间的发育特征，如孔隙尺寸、连通性、形态特征等9-11。孔隙网络模型自从 FATT 引入后，近年来有了长足的发展。赫瑞瓦特大学的研究团队12-15、TAH-MASEBI 团队16、孙建孟团队17、陶果团队18、刘向君团队19以及刘洁团队20都对建模方法进行了改进，并在各自的研究领域进行了推广。帝国理工大学BLUNT 团队7将微纳米计算机断层扫描应用到孔隙建模，使孔隙空间特性的直接测量得到了显著进步；林承焰等21在处理灰度图像时，采用了相比中值滤波法更为准确的非局部均值滤波算法(NLMfilter)进行滤波降噪处理；2 者使孔隙网络建模方法有了较大改进。笔者基于孔隙网络模型重构数字岩心，将综合孔隙网络建模方法结构性较强以及数字岩心建模方法所得数值精准的优点，研究不同图像处理以及孔隙网络提取方法对数字岩心重构结果的影响；之后将运用BLUNT 等22-23提出的孔隙级渗流数值模拟方法，研究孔隙网络模型提取方法对 CO2咸水层封存过程中多孔介质多相流过程的影响。1材料与方法1.1材料与设备1.1.1实验材料岩心选用储层较为常见的 Berea 砂岩，其直径为15.36mm，高度为40.77mm，渗透率为0.84041012m2；通过饱和称重法获得其孔隙度为 30.5%。1.1.2仪器设备采用型号为 InspeXioSMX-22CT 的微焦点 X 射线 CT 机作为获取 Berea 岩心 CT 图像的实验设备，其图像分辨率为 32m/pixel，最大视域为 300mm 直径和 300mm 高度，X 射线源的最大电压和电流为225kV 和 1mA11,24-25。该台式 CT 机主要包括：X 射线源、样品夹持器、X 射线探测器等部分，如图 1所示。ImageJ 软件是由美国国立卫生研究院开发的基于 Java 的图像处理程序；处理过的图像在 Matlab 平台进行数据的提取。等效孔隙网络可视化模型建立采用 Amira 软件实现。Amira 是澳大利亚 VisageImaging 公司研发的一款三维建模软件，用于数字岩心重构与可视化处理。1.2方法与过程1.2.1方法本次孔隙级渗流研究使用的方法主要有孔隙网络提取和两相渗流模型。(1)孔隙网络提取。最大球算法是由 SILIN 等在测量排水过程中无量纲毛细管压力时首先提出的26-27，之后 DONG等28-29改进其建模方法，提高了建模速度，促进了算法的成熟。近年来，最大球算法由于较高的建模准确性和较快的建模速度，得到了广泛应用30，因此笔者采用该方法提取了 Berea 岩心等效孔隙网络，主要步骤包括以下 3 步：搜索最大球。因为 3D 图像中的体素是离散和不连续的，很难给出精确的半径值，最大球算法中为了限定球的大小，设置了上限(Rright)和下限(Rleft)30。Rright和 Rleft的定义如式(1)和(2)，过程如图 2 所示。建立最大球连通关系。最大球算法的原理为：找出在图像中的每一个以体素为球心的最大内切球的球体，这个最大的球体只是与颗粒或边界接触，清除球体内的其他小球，这样描述出的网络空间没有冗余。从局部来看，最大球表示的是孔隙结构而大球与大球之间夹的小球表示的是喉道结构，其基本原理如图 3 所示11,30。孔隙与喉道的识别和参数计算。最大球算法通过设定一个球体半径比值(一般为 0.7)来划分孔隙体和喉道。如果一个最大球CT扫描仪控制电脑X射线投影机试料台图1Micro-CT 实验平台Fig.1Micro-CTexperimentalplatform(a)上限(b)下限图2最大球半径上下限示意28Fig.2Schematicdiagramofupperandlowerlimitsofmaximumsphericalradius282804煤炭学报2023年第48卷的半径是其上一代最大球的 0.7 倍及以上时，被识别为孔隙体，反之，则被认为是喉道30。R2right=dist2(C0,Vg)=(xgxc)2+(ygyc)2+(zgzc)2,Vg Sg,C0 S0(1)R2left=maxdist2(V0,C0)?dist2(V0,C0)R2right,V0 S0,C0 S0(2)式中，Vg为最近颗粒体素，(xg，yg，zg)为其空间坐标；V0为孔隙体素；C0为孔隙中心素，(xc，yc，zc)为其空间坐标；S0和 Sg分别为孔隙和固体体素点。孔隙网络的各项参数借助 DONG和BLUNT 开发的程序28-29来计算。其中，孔隙或喉道的体积大小取决于每个块的孔隙或喉占用的体素的数量；配位数由与孔隙体相连的喉道数目来决定；孔隙和喉道结构的不规则程度借用无量纲参数形状因子 G 来表征11，其物理定义为孔隙和喉道结构的横截面积 As与周长Ps平方的比值，如式(3)所示。根据形状因子 G，可以用具有等同几何特征的简单几何体如三角形、正方形以及圆形来代替实际孔隙或喉道，如图 4 所示11。G=AsP2s(3)图 2 演示了定义最大球半径的一个简单过程，Rright定义的球体半径是 8 个小球(02+22+22=8)，而Rleft定义的球体半径是 6 个小球(12+12+22=6)。其定义的最大球体半径 R0取值范围是(Rleft，Rright)11。(2)两相渗流模型。帝国理工大学 VALVATNEL31开发的孔隙网络两相渗流模拟程序可以分析预测流动的相对渗透率和毛管压力，被广为应用，下面是模型的基本情况。模型的基本假设。CO2和咸水均为不可压缩的牛顿流体，且渗流过程中 2 者互不相溶；忽略黏性力和毛细管压降的影响，流动主要由毛细管力控制25。模型对润湿滞后效应的处理。对于润湿滞后效应的处理采用 MORROW32通过大量的实验得到的润湿滞后模型，如图 5 所示25。模型对润湿性的处理。上述模型可以模拟单一润湿和混合润湿 2 种润湿情况，笔者对于 CO2咸水的模拟使用单一润湿系统。在单一润湿系统中气水接触角的变化会改变润湿特性，具体表现为：当接触角增大时，Berea 岩心的润湿特性由原来的中性亲水逐渐开始变得亲气，如图 6 所示25。孔隙级渗流模拟的基本过程包括驱替过程和吸球AR42RRR3R1R2家庭A家庭B共同的孩子祖宗A祖宗B图3最大球算法示意28Fig.3Schematicdiagramofmaximumballmethod28G=AsPs2116G=G=(0,G=14336rrrrAsPs图4截面几何形状随形状因子变化Fig.4Graphofgeometricshapechangewithshapefactor0306090120150180306090120150180后退接触角前进接触角本征接触角/()后退/前进触角/()图5润湿滞后模型示意Fig.5Schematicdiagramofwettinghysteresismodel第7期美合日阿依穆太力普等：基于孔隙网络模型的气液渗流特性2805入过程。借助该模型，驱替、吸入过程中均可以计算得到气水两相的相对渗透率和毛细管压力随着饱和度的变化11。1)驱替过程。为了模拟结果的准确性，在驱替开始前需要把已经用水饱和的孔隙网络模型的入口端与气源相连接，因此整个过程中水相的压力保持不变。开始驱替后，随着气体的进入，气相压力逐步升高，气相将把该单元中的水驱走。当模型中含水饱和度达到束缚水饱和度时，该驱替过程结束11。2)吸入过程。吸入过程中指的是液相水重新占领气体占据的孔隙空间，完成水对气的驱替，主要包括活塞式驱替、孔隙体充填和卡断或阶跃 3 种方式25。吸入过程结束后，再次把孔隙网络模型的入口端与驱替相流体(气)源相连接，可以连续进行二次驱替模拟。因此，与物理实验相比，孔隙级渗流模拟可以多次模拟渗流过程，是进行渗流理论研究的一种快速、经济的技术手段25。3)参数计算。饱和度计算。通过计算出每一孔隙、喉道中的气水量后，统计所有孔隙、喉道中的含水量而求得整个孔隙网络模型的含水饱和度 Sw：Sw=ni=1Viw/ni=1Vi(4)式中，n 为孔隙网络模型中所有喉道和孔隙的总数；Vi为第 i 个喉道(或孔隙)的体积；Viw为第 i 个喉道(或孔隙)中水所占的体积。毛细管压力计算。(a)驱替过程。考虑两相界面间力的平衡原则，驱替过程中毛细管入口压力 pC计算如下：pC=cosr(1+2G)rFd(r,G,)(5)式中，为气液界面张力；Fd为无量纲修正因子，它是后退角 r、形状因子 G 和多边形内半角 的函数；r为内切半径。(b)吸入过程。在自发侵入过程中，孔隙体的毛细管进入压力将取决于相邻充油喉道的数量 n(协同孔隙体填充)。由于充油喉道的准确空间位置难以估计，毛细管入口压力通常采用式(6)进行计算：pC=2cosarni=1Aixi(6)A2An=0.03/k式中，a为前进角，()；n 为与孔隙相连的喉道中充满驱体相油/气的喉道数；Ai为任意参数，其与绝对渗透率 k 的关系为，因此单位为 m1，且定义 A1=0；xi为 01 之间的随机数。渗透率的计算。模型的绝对渗透率是基于达西公式得到的。其过程为：先将整个模型饱和一种流体，之后给模型施加一个驱动压力(PinPout)，最后统计流体流量算出模型的绝对渗透率：k=pqispLA(PinPout)(7)式中，k 为绝对渗透率，m2；p为 p 相流体的黏度，Pas；qisp为孔隙网络模型完全饱和 p 相流体时的总流量，m3/s；L 为模型长度，m；A 为模型面积，m2。当模型饱和两相流体时，通过某相流体在多相流时的流量与单相时的流量之比可以计算该相流体的相对渗透率 krp，即krp=qtmpqtsp(8)式中，qtmp为 p 相在多相流时的流量；qtsp为单相流时的流量。1.2.2建模过程先是基于 micro-CT 技术，扫描得到 Berea 岩心的 CT 图像；然后采用 ImageJ 软件，通过尺寸调整、滤波降噪、阈值分割等一系列图像处理过程，得到了Berea 岩心的黑白图像；之后对应函数计算将处理过的图像在 MATLAB 平台转化为二值化数据；然后基于最大球算法将孔隙空间图像转化成孔隙和喉道的简化等效网络；最后采用 Amira 软件建立能够生动描述 Berea 岩心孔隙网络的球杆模型，从而得到孔隙网络模型，具体过程如图 7 所示。2参数统计与分析2.1建模参数笔者使用 Berea 岩心作为实验样品，由 CT 实验得到了 Berea 岩心 CT 图像，如图 8 所示。首先，裁剪获得 Berea 岩心的 CT 图像的有效区域。其次，通过Mean 和 Median 两种滤波方法对图像进行处理，分别水气亲水孔隙壁面亲气孔隙壁面气体运移图6润湿性转变示意Fig.6Schematicdiagramofwettabilityalternation2806煤炭学报2023年第48卷将其半径设置为 r1、r2。之后，分别采用 IsoData、Mean、Otsu 以及 Shanbhag 四种阈值分割的算法对图像进行处理。最后，对黑白图像进行二值化处理，并基于最大球算法提取岩心等效孔隙网络，采用 Amira软件实现球杆网络模型的可视化，重构 Berea 岩心可视化孔隙网络模型。笔者基于岩心 micro-CT 图像的数字岩心重构所选取的相关算法与参数见表 1。2.2模型特征参数分析基于不同图像处理算法以及参数得到了 Berea 岩心等效孔隙网络孔隙度、连通性(配位数)、几何尺寸(孔隙和吼道的半径与体积)和形状因子等特征参数，并对其进行统计与对比分析，揭示不同图像处理算法以及参数对基于最大球算法提取岩心等效孔隙网络结果的影响机制。表 2 列出了不同方法得到的网络模型各项参数，包括孔隙度、孔隙数目、喉道数目、长度、体积等特征参数。首先，对模型的孔隙度进行比较，对比结果如图 9(a)所示。与其他模型相比，Pic1、Pic2、Pic3 处理方式与实际 Berea 岩心孔隙度实验值偏差较大，其中Pic2 和 Pic3 的偏差高达 32.47%。之后，笔者选取孔隙网络提取方法更具代表性且结果较可靠的 Pic4、Pic5、Pic6、Pic7 四个模型，对表征其连通性的配位数，表示孔隙空间内粒度分布的孔隙与喉道几何参数包括了半径、体积和长度，以及其形状因子等特征参数进行统计和对比分析。2.2.1配位数配位数(连接数)是指与孔隙相连的独立的喉道数，能够度量孔隙网络连通性，是评估多孔介质运输性质重要特征参数33。对不同方法得到的配位数进行对比分析，结果如图 9(b)所示。由图 9(b)可知，各模型的配位数分布在 30 以内且均呈负指数分布，峰值在 5 左右，这是由于这 4 种模型均来自同一 Berea二值化数据球杆模型Berea岩石样品CT图像二值化图像孔隙网络模型参数Micro-CTB1B2S1S2L1获取图像图像处理Image JavaAmira重建模型最大球算法提取孔隙网络Matlab数据提取喉道半径/m喉道形状因子 喉道总长度/m1.689 61051.143 61056.613 51052.880 01051.301 01058.446 91051.902 31052.900 01050.022 7990.026 6290.020 3440.031 2310.030 0500.032 8110.027 9070.037 7968.683 01062.029 31048.683 01068.769 41053.908 41048.683 01067.910 61051.833 9104图7孔隙网络模型的构建流程Fig.7Constructionprocessofporenetworkmodel(a)截面原始 CT 图像(b)正向原始 CT 图像图8Berea 岩心 CT 图像Fig.8CTimageofBereacore表 1 基于 8 种不同方法对 Berea 岩心 CT 图像进行处理Table 1 Image processing of Berea sandstone CT based oneight different methods模型图像处理模型提取滤波半径/pixels分割阈值最小球半径/voxelsPic1r1=2,r2=1.5IsoData2Pic2r1=2,r2=1.5Mean1Pic3r1=2,r2=1.5Mean2Pic4r1=1.5,r2=2IsoData2Pic5r1=1.5,r2=2Shanbhag2Pic6r1=2,r2=1.5Shanbhag2Pic7r1=1.5,r2=2Otsu1Pic8r1=1.5,r2=2Otsu2第7期美合日阿依穆太力普等：基于孔隙网络模型的气液渗流特性2807岩心样本，有着相同的粒径分布形态。然而，Pic4 和Pic5 模型在峰值处频率较低，分布较为分散，Pic6 的分布最为集中，说明不同的阈值分割方式对配位数有较为明显的影响。2.2.2孔喉几何尺寸的大小孔隙和喉道的尺寸主要是指孔隙和喉道的半径。从图 10 可以看出，总体上孔隙半径分布在 80m 内，呈中间高两边低的分布态势，孔隙半径集中在 20m左右；吼道半径分布在 60m 内，呈负指数递减的分布态势，喉道半径集中在 10m 左右，说明其均质性较好。其中 Pic6 和 Pic7 峰值更高，较其他方式分布更集中。由图 11 可知，4 种方法提取的孔隙体积和喉道体积分布情况基本一致，均呈负指数递减的分布态势。孔隙体积分布在 1.251011m3以内，峰值在 0.51011m3；吼道体积分布在 81013m3以内，峰值在21013m3。而喉道体积 Pic4 在峰值处最小，说明其得到的吼道体积偏大，这可能是由于 isodata 这种阈值分割方法识别出的孔喉偏大导致的。图 12 分别为喉道的总长度 lij和喉道的有效长度lt的分析对比。从图 12 可以看出，喉道的总长度 lij和喉道的有效长度 lt呈负指数递减的分布态势。喉道的总长度 lij分布在 1mm 以内，峰值在 0.2mm；喉道的有效长度 lt分布在 0.15mm 以内，峰值在 0.03mm。表 2 Berea 岩心提取的孔隙网络特征参数Table 2 Extracted pore network properties of Berea sandstone模型Pic1Pic2Pic3Pic4Pic5Pic6Pic7Pic8孔隙度/%24.87440.40340.40327.48427.96830.26627.00327.003孔隙数量/个1552133412219911498121417661515喉道数量/个26573519339518062805242128742686平均配位数3.305.135.403.503.603.803.103.4最大配位数1728281819182120孔隙平均半径/m23.026.028.028.724.027.321.423.6孔隙最大半径/m60.769.871.272.572.565.967.973.3孔隙平均体积/1012m30.51.11.21.00.70.90.60.6孔隙最大体积/1012m36.814.113.910.97.09.37.57.3孔隙平均形状因子/1023.12.82.73.02.92.93.02.9喉道平均半径/m10.212.612.713.510.812.310.310.6喉道最大半径/m49.552.250.861.654.559.659.752.4喉道平均体积/1015m336.232.634.160.035.945.836.039.3喉道最大体积/1015m3513.2803.3796.7985.9733.22775.11639.21637.9喉道平均总长度/m170.3219.1223.1208.4179.6195.8173.9178.9喉道平均有效长度/m25.931.432.032.227.029.226.226.9喉道平均形状因子/1023.153.153.133.123.143.143.143.10010203040Pic6Pic8Pic7Pic5Pic4Pic3Pic2计算值Pic1模型实验值1.00.53020100配位数(a)(b)Pic4 Pic5 Pic6 Pic7孔隙度/%归一化频率图9不同模型孔隙度及配位数分布对比Fig.9Comparisonofporosityofandcoordinationnumberdistributionofdifferentsamples2808煤炭学报2023年第48卷模型 Pic4 的喉道的总长度 lij分布峰值不明显，数值不集中；模型 Pic7 和 Pic6 的提取结果基本相似。由此可知，阈值分割采用 Shanbhag 算法或 Otsu 算法所提取的孔隙网络喉道长度基本相同。2.2.3形状因子以孔隙网络模型为基础进行孔隙级渗流数值模拟时，形状因子是很重要的特征参数。图 13 表示对数字岩心孔隙与喉道的形状因子进行统计和对比分析。总体上 4 种方法得到的孔隙形状因子和吼道形状因子分布形态相当，呈单峰分布的态势，分布在0.07 以内。孔隙形状因子峰值在 0.03 左右，吼道形状因子峰值在 0.035 左右，进而可以看出该岩心实际孔隙与吼道形状近似三角形，与 Berea 岩心实际情况较为吻合。这是由于所有方法是基于同一个 Berea 岩心样品进行的进一步处理，因此它们的颗粒排布一致，形状因子的分布也因此几乎一致。由图 13 可知，Pic5 得到的孔隙形状因子最集中，数值偏小，而 Pic6得到的形状因子峰值最低，数值偏大，这可能是不同的滤波处理方式造成的。3孔隙级渗流模拟为了进一步研究 Berea 岩心孔隙网络模型的提取(a)(b)0.30.40.20.100.080.060.04孔隙半径/mm喉道半径/mm Pic4 Pic5 Pic6 Pic7 Pic4 Pic5 Pic6 Pic70.020.60.50.40.30.20.100.060.040.02归一化频率归一化频率图104 个模型孔隙及喉道的半径分布对比Fig.10Comparisonofdistributionofporeandthroatradiusoffoursamples(a)(b)孔隙体积/1012 m3喉道体积/1013 m31.00.80.60.40.20151051.00.80.60.40.2015105归一化频率归一化频率 Pic4 Pic5 Pic6 Pic7 Pic4 Pic5 Pic6 Pic7图114 个模型孔隙及喉道的体积分布对比Fig.11Comparisonofporeandthroatvolumedistributionoffoursamples(a)(b)喉道总长度/mm喉道有效长度/mm Pic4 Pic5 Pic6 Pic7 Pic4 Pic5 Pic6 Pic70.80.60.40.200.80.60.4归一化频率0.20.80.60.40.200.120.150.090.06归一化频率0.03图124 个模型的喉道总长度和有效长度的分布对比Fig.12Comparisonofthedistributionoftotalthroatlengthandeffectivethroatlengthoffoursamples第7期美合日阿依穆太力普等：基于孔隙网络模型的气液渗流特性2809过程与计算方法对 CO2咸水(0.103mol/kgNaCl 溶液)层封存过程中气液在多孔介质中的多相流过程，笔者基于孔隙级数值模拟方法预测了 Pic4 和 Pic7 等效孔隙网络中两相流动过程，并与实验值进行了对比。孔隙级渗流数值模拟主要输入参数见表 3。CO2驱替咸水过程中，气液两相相对渗透率随含水饱和度的变化过程如图 14(a)、(c)所示，从图中可知，Pic4 和 Pic7 的预测曲线有较大差异，Pic7 的 CO2相对渗透率较 Pic4 的大，而咸水相对渗透率较 Pic4的小，吸入过程也是如此。这是由于 Pic7 的孔喉半径比 Pic4 的小，通过 Pic7 的方法提取出岩心更小的孔隙和喉道。且 Pic7 的咸水相对渗透率与实验值较接近，而 Pic4 的 CO2相对渗透率与实验值较接近。总的来说，在驱替和吸入过程中 Pic7 的毛细管压力比(a)(b)归一化频率归一化频率0.60.30.060孔隙形状因子0.030.30.20.10.070.060.050.040.030