本发明公开一种基于孔隙度重构数字岩心的原位温压两相流分析方法。采用恒温岩心夹持器提供地层原位温压；先后获取流体F1饱和岩心CT图像IM1、流体F2饱和岩心CT图像IM2；通过IM1和IM2图像配准，计算差值图像IM21，通过图像分割提取流体相，计算两种流体灰度差μF21；用IM21除μF21获得孔隙度图像IMε，利用孔渗关系模型，获得渗透率图像IMκ；利用流体热力学软件，获取不同温压下流体F3和F4的粘度和密度函数；将岩心孔隙度和渗透率图像、流体粘度和密度函数导入CFD软件，开展两相流仿真，获得流体流速、压力和含量分布。以孔隙度重构替代孔裂隙空间重构，降低成像分辨率要求，为原位温压下岩心两相流提供可视化研究，实施步骤简单、成本低。
1.一种基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于步骤如下：
步骤1. 将设有导水钻孔的岩心放入恒温岩心夹持器，采用外部缠绕加热恒温套的岩心夹持器模拟岩心所在地层的原位温度、围压、轴压和孔压环境；
步骤2. 选择X射线衰减系数差异大的流体F1和流体F2先后饱和岩心，利用扇束X射线成像系统使用在相同X射线源能量和相近的成像分辨率条件下，获取流体F1饱和岩心CT图像IM1、流体F2饱和岩心CT图像IM2；
步骤3. 选择以流体F1饱和岩心CT图像IM1为模板，对流体F2饱和岩心CT图像IM2进行图像配准，获得配准后的图像IM2’，通过用图像IM2’减去图像IM1计算体现岩心中流体变化的差值图像IM21，通过图像分割提取流体F1和流体F2，计算两种流体的平均灰度差μF21；
步骤4. 用差值图像IM21除以流体平均灰度差μF21获得整个岩心的孔隙度图像IMε，利用Kozeny‑Carman孔渗关系模型，由孔隙度图像经像素运算获得整个岩心的渗透率图像IMκ；
步骤5. 利用流体热力学软件获取不同温压下流体F3和流体F4以温度和压力为自变量的粘度和密度函数；
步骤6. 将岩心ROI区域的孔隙度图像和渗透率图像、流体粘度和密度函数导入CFD软件，开展流体F3和流体F4的两相流仿真，获得流体F3和流体F4在岩心ROI中的渗流流速、压力和流体含量分布。
2.根据权利要求1所述的基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：所述的恒温岩心夹持器，包括三轴岩心夹持器，三轴岩心夹持器内设有岩心/围压腔、轴压腔、孔压入口及出口，三轴岩心夹持器外侧设有柔性加热套，三轴岩心夹持器的围压腔、轴压腔、孔压入口及出口分别设有围压阀、轴压阀、入口孔压阀及出口孔压阀，夹持器外壳为低X射线衰减材料，对于热变形可忽略的岩心不用加热恒温套。
3.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：流体F1和流体F2，须选择不与岩心反应、且X射线衰减系数差异大的流体，同时流体F1应为低残留流体。
4.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：所述岩心的一个端部设有导水钻孔，导水钻孔为大于Ф10像素×10像素的圆柱形导水钻孔。
5.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：采用阈值算法进行图像分割，从图像IM1和图像IM2’图像中的导水钻孔处分离流体F1和流体F2，采用像素平均算法，计算流体F1的灰度平均值μF1以及流体F2的灰度平均值μF2，由μF2减去μF1得到流体平均灰度差μF21。
6.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：基于Kozeny‑Carman孔渗关系模型或者单相渗流实验拟合得出的渗透率‑孔隙度关系模型κ = f(ε)，对孔隙度图像IMε开展像素运算获得渗透率图像IMκ，即IMκ(i,j,k) = f(IMε(i,j,k))，Kozeny‑Carman孔渗关系模型如下：
式中，κ—渗透率；ε—孔隙度；κR—单相渗流实验获得的岩心渗透率；εR—压汞实验获得的岩心孔隙度。
7.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：所述的岩心ROI区域的获取方法为，在IMε或IMκ图像上建立与岩心同轴且略小于岩心的圆柱选区，将选区外部的像素赋值为0，然后通过图像裁剪将去除人工导水孔段的岩心。
8.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：所述步骤6中的孔隙度和渗透率图像导入CFD软件的方法为，采用插值算法，由IMε(i,j,k)和IMκ(i,j,k)的像素点集分别建立孔隙度和渗透率的空间分布ε(x,y,z)和κ(x,y,z)。
9.根据权利要求1所述基于孔隙度重构数字岩心的原位温压两相流分析方法，其特征在于：所述步骤6中的两相流仿真，其具体步骤如下：基于孔隙度和渗透率的空间分布建立原位温压下岩心的几何和材料模型，基于流体F3和流体F4的粘度和密度函数建立原位温压下两种流体的材料模型，采用六面体或四面体单元对几何模型进行网格划分，以Cahn–Hilliard–Brinkman耦合方程组作为控制方程，以岩心原位孔隙水压为出口压力边界，以原位流流速或体流量为入口边界，求解单元节点及单元平均流速、压力和流体含量，Cahn–Hilliard–Brinkman耦合方程组作为控制方程如下：
式中，ρ1—流体F3平均密度，随孔隙水压和温度变化，ρ2—流体F4平均密度，随孔隙水压和温度变化；ε—孔隙度；u—渗流流速；t—时间；p—压力；I—单位矩阵；μ1—流体F3平均粘度，随孔隙水压和温度变化，μ2—流体F4平均粘度，随孔隙水压和温度变化；κ—渗透率；β—Forchheimer系数；σ—界面张力；ϕ—相变量；α—界面厚度控制量；χ—流动性调定参数；
ψ—相场方程辅助变量；下标1、2代表流体F3和流体F4两种流体。
一种基于孔隙度重构数字岩心的原位温压两相流分析方法
技术领域
[0001] 本发明属于二氧化碳地质封存与利用、数字岩心领域，具体涉及一种基于孔隙度重构数字岩心的原位温压两相流分析方法。
背景技术
[0002] 二氧化碳地质封存与利用技术，是指将大型排放源捕集的CO2注入咸水层、油气田、煤层、盐床、地热田等深部地层，进行长期封存，以实现CO2减排和地下资源增产的过程。二氧化碳储层深度一般大于800米，地温大于31.1℃，地应力大于7.38MPa。深部岩石通常具有孔隙尺寸小、渗透率低、非均质性强的特点。而且，深地高温高压环境下岩心的孔裂隙结构与没有原位温压约束的岩心孔裂隙结构存在差异，深地原位温压下流体的特性与常温常压下流体的特性可能存在巨大差异，尤其是CO2。高温、高孔隙水压和高地应力环境下非均质致密岩心中的两相流(CO2－咸水、CO2－油、CO2－甲烷等)规律是二氧化碳地质封存与利用技术中的一个热点课题。
[0003] 基于孔隙空间重构的数字岩心方法，已成功应用于多孔介质渗流的研究中。高分辨X射线成像被用于获取精确的孔裂隙结构，通过二值化建立三维孔裂隙空间，通过Navier‑Stokes模型模拟流体在三维孔隙空间中的流动，获得三维孔裂隙空间的渗流特性。然而该方法中岩心X射线成像通常在常温常压下进行。随着4D CT技术的发展，通过高分辨显微成像结合高温高压岩心渗流实验，可在线捕获流体在孔隙空间内流动的图像，直接用于分析渗流特性。但是，4D CT实时成像存在设施少、成本高、原位温压渗流实验架设难度大的问题。此外，两种方法都无法避开分辨率成像与样品尺寸的矛盾，高分辨率意味着大尺寸孔隙或者小尺寸样品。例如，平均孔径50μm的均质砂岩，若要获取精确的孔隙结构，其成像分辨率最大不宜超过10μm，则岩心直径最大不能超过1.024cm；对于致密岩心，则需要更小的直径；对于非均质致密岩心，高分辨X射线成像可接受的样品尺寸通常不具代表性。
发明内容
[0004] 针对现有技术的不足之处，提供一种实施步骤简单、成本低的一种基于孔隙度重构的原位温压下岩心两相流分析方法。
[0005] 为实现上述技术目的，本发明的一种基于孔隙度重构的原位温压下岩心两相流分析方法，步骤如下：
[0006] 步骤1.将设有导水钻孔的岩心放入恒温岩心夹持器，采用外部缠绕加热恒温套的岩心夹持器模拟岩心所在地层的原位温度、围压、轴压和孔压环境；
[0007] 步骤2.选择X射线衰减系数差异大的流体F1和流体F2先后饱和岩心，利用扇束X射线成像系统使用在相同X射线源能量和相近的成像分辨率条件下，获取流体F1饱和岩心CT图像IM1、流体F2饱和岩心CT图像IM2；
[0008] 步骤3.选择以流体F1饱和岩心CT图像IM1为模板，对流体F2和岩心CT图像IM2进行图像配准，获得配准后的图像IM2’，通过用图像IM2’减去图像IM1计算体现岩心中流体变化的差值图像IM21，通过图像分割提取流体F1和流体F2，计算两种流体的平均灰度差μF21；
[0009] 步骤4.用差值图像IM21除以流体平均灰度差μ21获得整个岩心的孔隙度图像IMε，利用Kozeny‑Carman孔渗关系模型，由孔隙度图像经像素运算获得整个岩心的渗透率图像IMκ；
[0010] 步骤5.利用流体热力学软件获取不同温压下流体F3和流体F4以温度和压力为自变量的粘度和密度函数；
[0011] 步骤6.将岩心ROI区域的孔隙度图像和渗透率图像、流体粘度和密度函数导入CFD软件，开展流体F3和流体F4的两相流仿真，获得流体F3和流体F4在岩心ROI中的渗流流速、压力和流体含量分布。
[0012] 进一步，所述的恒温岩心夹持器，包括三轴岩心夹持器，三轴岩心夹持器内设有岩心/围压腔、轴压腔、孔压入口及出口，三轴岩心夹持器外侧设有柔性加热套，三轴岩心夹持器的围压腔、轴压腔、孔压入口及出口分别设有围压阀、轴压阀、入口孔压阀及出口孔压阀，夹持器外壳为低X射线衰减材料，对于热变形可忽略的岩心可以不用加热恒温套。
[0013] 进一步，流体F1和流体F2，须选择不与岩心反应、且X射线衰减系数差异大的流体，同时流体F1应为低残留流体，优选地，流体F1为空气、He或N2气体，流体F2为水。
[0014] 进一步，所述岩心的一个端部设有导水钻孔，导水钻孔为大于Ф10像素×10像素的圆柱形导水钻孔。
[0015] 进一步，采用阈值算法进行图像分割，从图像IM1和图像IM2’图像中的导水钻孔处分离流体F1和流体F2，采用像素平均算法，计算流体F1的灰度平均值μF1以及流体F2的灰度平均值μF2，由μF2减去μF1得到流体平均灰度差μF21。
[0016] 进一步，基于Kozeny‑Carman孔渗关系模型或者单相渗流实验拟合得出的渗透率‑孔隙度关系模型κ＝f(ε)，对孔隙度图像IMε开展像素运算获得渗透率图像IMκ，即IMκ(i,j,k)＝f(IMε(i,j,k))，
[0017] Kozeny‑Carman孔渗关系模型如下：
[0018]
[0019] 式中，κ—渗透率；ε—孔隙度；κR—单相渗流实验获得的岩心渗透率；εR—压汞实验获得的岩心孔隙度。
[0020] 进一步，所述的岩心ROI区域的获取方法为，在IMε或IMκ图像上建立与岩心同轴且略小于岩心的圆柱选区，将选区外部的像素赋值为0，然后通过图像裁剪将去除人工导水孔段的岩心。
[0021] 进一步，所述步骤3‑4中对图像配准、差值、图像分割、像素运算，使用的图像处理软件包括IMAGEJ、MATLAB Image Processing Toolbox、AVIZO、SIMPLEWARE、DRAGONFLY；
[0022] 所述步骤5中的流体热力学软件包括COCO/TEA、COOLPROP、REFPROP；
[0023] 所述步骤6中的CFD软件可采用FLUENT、COMSOL、OPENFORM；
[0024] 所述步骤6中的孔隙度和渗透率图像导入CFD软件的方法为，采用插值算法，由IMε(i,j,k)和IMκ(i,j,k)的像素点集分别建立孔隙度和渗透率的空间分布ε(x,y,z)和κ(x,y,z)。
[0025] 进一步，所述步骤6中的两相流仿真，其具体步骤如下：
[0026] 基于孔隙度和渗透率的空间分布建立原位温压下岩心的几何和材料模型，基于流体F3和流体F4的粘度和密度函数建立原位温压下两种流体的材料模型，采用六面体或四面体单元对几何模型进行网格划分，以Cahn–Hilliard–Brinkman耦合方程组作为控制方程，以岩心原位孔隙水压为出口压力边界，以原位流流速或体流量为入口边界，求解单元节点及单元平均流速、压力和流体含量，
[0027] Cahn–Hilliard–Brinkman耦合方程组作为控制方程如下：
[0028]
[0029] 式中，ρ1—流体F3平均密度，随孔隙水压和温度变化，ρ2—流体F4平均密度，随孔隙水压和温度变化；ε—孔隙度；u—渗流流速；t—时间；p—压力；I—单位矩阵；μ1—流体F3平均粘度，随孔隙水压和温度变化，μ2—流体F4平均粘度，随孔隙水压和温度变化；κ—渗透率；β—Forchheimer系数；σ—界面张力；φ—相变量；α—界面厚度控制量；χ—流动性调定参数；ψ—相场方程辅助变量；下标1、2代表流体F3和流体F4两种流体。
[0030] 有益效果：
[0031] 1)以往基于孔裂隙空间重构的数字岩心分析方法，需要通过高分辨X射线成像，精确解析孔裂隙空间，大大限制了样品尺寸；尤其对于深地非均质致密岩心，由于岩心孔径较小，只能使用微型样品，从而造成分析结果代表性差；本发明基于孔隙度重构的数字岩心方法——即利用两种X射线衰减系数差异较大的流体去饱和同一块岩心，通过两种流体饱和岩心的X射线CT图像的差异，重构岩心孔隙度图像——不需要精确解析孔隙空间，大大降低对成像分辨率的要求，解除了样品尺寸的限制，适用于大部分实验室尺度岩心及标准岩心。
[0032] 2)地层原位温压下岩心两相流的可视化研究主要依赖4D CT技术，然而目前4D CT实时成像设施少、成本高、原位温压渗流实验架设难度大，同样依赖高分辨X射线成像以获得精确解析的流体界面和孔隙空间；本发明通过原位温压下岩心X射线成像、岩心孔隙度重构和数值仿真相结合，不仅对成像分辨率要求低，更不需要在成像过程中开展渗流实验，从而大大降低了岩心原位温压两相流可视化研究的实验难度和成本。
附图说明
[0033] 图1为本发明基于孔隙度重构数字岩心的原位温压两相流分析方法流程示意图；
[0034] 图2为本发明模拟原位温度、围压、轴压和孔压环境的恒温岩心夹持器结构示意图；
[0035] 图3为本发明使用的扇束X射线成像系统示意图；
[0036] 图4为本发明恒温岩心夹持器内岩心CT正交视图；
[0037] 图5为本发明中整个岩心孔隙度重构图像(IM2’‑IM1)/(μ2‑μ1)示意图；
[0038] 图6为本发明中不同温压下水和CO2流体的粘度和密度函数示意图；
[0039] 图7为孔隙度和渗透率图像导入CFD软件中建立的岩心几何、网格及材料模型示意图；
[0040] 图8为CO2驱水过程中CO2和水的压力、流速及含量空间分布的模拟结果示意图。
具体实施方式
[0041] 下面结合附图对本发明的实施例作进一步说明：
[0042] 实施例一、
[0043] 一种基于孔隙度重构数字岩心的原位温压两相流分析方法，具体步骤如下：
[0044] 步骤1.采用如图2所示恒温岩心夹持器，将Ф50×160mm的砂岩岩心装入夹持器中，岩心一端导水钻孔尺寸为Ф7.5×10mm，打开轴压阀和围压阀，利用恒压泵以液压油为介质施加5MPa轴压和17MPa围压环境，孔压出口阀常闭，打开孔压入口阀，以空气为介质施加14MPa孔压环境，通过柔性加热保温套提供65℃温度环境，待温压稳定48h后认为岩心饱和，关闭轴压阀、围压阀以及孔压入口阀，恒温岩心夹持器脱离恒压泵；
[0045] 步骤2.采用图3所示扇束X射线成像系统，在135KeV/200mA的射线能量和0.3mm/像素分辨率条件下获取空气饱和岩心CT图像IM1，然后换水作为孔压介质，在相同温压环境下饱和相同时间后，用同一套扇束X射线成像系统，在相同射线能量和分辨率条件下，获取水饱和岩心CT图像IM2；
[0046] 步骤3.以空气饱和岩心CT图像IM1为模板，对水饱和岩心CT图像IM2进行图像配准，采用ImageJ图像处理软件进行图像配准，获得配准后的图像IM2’，原始IM1、IM2图像和配准后IM2’图像的正交视图如图4所示，图4中分别体现空气饱和岩心CT正交视图IM1，水饱和岩心CT正交视图IM2，配准后的水饱和岩心CT正交视图IM2’，图4中A1和B1分别表示空气饱和岩心的两端CT视图，A2和B2分别表示水饱和岩心的两端CT视图，A3和B3分别表示配准后的水饱和岩心的两端CT视图，序号1表示岩心，序号2表示导水钻孔处空气，序号3表示导水钻孔处水，序号4表示胶套管序号，5表示围压舱液压油，序号6表示夹持器内衬，序号7表示恒温岩心夹持器的夹持器外壳；
[0047] 用正交视图IM2’减去正交视图IM1，获得差值图像IM21，采用ImageJ进行阈值分割，从空气饱和岩心IM1分离出导水钻孔处的空气(图4中序号2)，得到空气的平均灰度值，从水饱和岩心CT图像IM2’分离出导水钻孔处的水(图4中序号3)，计算水的平均灰度值，以空气的平均灰度值减去水的平均灰度值，获得流体平均灰度差μ21＝725；
[0048] 步骤4.采用ImageJ，以图像IM21除以流体平均灰度差μ21，获得岩心孔隙度图像IMε，孔隙度值在0‑0.4之间，如图5所示，利用Kozeny‑Carman孔渗关系模型，由孔隙度图像经像素运算获得渗透率图像IMκ；即IMκ(i,j,k)＝f(IMε(i,j,k))，
[0049] Kozeny‑Carman孔渗关系模型如下：
[0050]
[0051] 式中，κ—渗透率；ε—孔隙度；κR—单相渗流实验获得的岩心渗透率；εR—压汞实验获得的岩心孔隙度；
[0052] 步骤5.利用COOLPROP开展流体热力学仿真，得到不同温度T和压力p下CO2和水的粘度和密度函数，以温度和压力为自变量，绘制CO2和水的粘度和密度分布如图6所示；
[0053] 步骤6.如图5所示，采用线性差值，将岩心ROI区域的孔隙度和渗透率图像导入Comsol Multiphysic软件建立岩心几何和材料模型，采用六面体单元对几何模型进行网格划分，岩心几何、网格及材料模型如图7所示，将流体粘度和密度函数代入Comsol Multiphysic软件建立CO2和水在65℃环境下随孔压变化的材料模型，以Cahn–Hilliard–Brinkman耦合方程组作为控制方程，设置岩心一端为入口流量边界，流量值为1.72ml/min，另一端为出口压力边界，出口背压值为14MPa，开展两相流仿真，获得CO2和水的流速、压力和含量分布图，如图8所示。
[0054] Cahn–Hilliard–Brinkman耦合方程组作为控制方程如下：
[0055]
[0056] 式中，ρ1—流体F3平均密度，随孔隙水压和温度变化，ρ2—流体F4平均密度，随孔隙水压和温度变化；ε—孔隙度；u—渗流流速；t—时间；p—压力；I—单位矩阵；μ1—流体F3平均粘度，随孔隙水压和温度变化，μ2—流体F4平均粘度，随孔隙水压和温度变化；κ—渗透率；β—Forchheimer系数；σ—界面张力；φ—相变量；α—界面厚度控制量；χ—流动性调定参数；ψ—相场方程辅助变量；下标1、2代表流体F3和流体F4两种流体。
