高密度电性参数监测的多尺度岩芯渗透率测定技术的制作方法

岩芯渗透率是描述岩石让液体透过能力大小的重要参数，它通过对岩芯的测定而获得，是水文地质、地质工程、石油地质、水文水资源、土壤学、矿产资源等领域的重要基础性参数，是后续模拟计算的重要输入变量。在室内和原位有效测定岩芯渗透率不仅一直是科学研究的一项重要内容，同时也是各种工程应用的迫切需要。本发明提出了一种基于高密度的电性参数监测获得不同尺度和方向上的岩芯渗透率参数的技术，在相关领域的科学研究和工程实践中都具有广泛的应用意义。

背景技术：

岩芯渗透率一般在毫达西量级，对其测定目前主要有定水位法和变水位法以及定流量法和压力脉冲法。前者通过在试件的两端施加一定的或变化的水压差，通过测量渗透流量来计算试件的渗透率参数。后者则在试件的一端以一定的流量注水或直接施加压力脉冲，通过测量试件两端间压力差随时间的变化来计算试件的渗透率参数。但无论是前者还是后者，都以测定试件两端的压力（或压力差）为主，再通过一维的流体动力学方程的求解及其和实验结果的拟合来获得参数。因此，所得参数是一个关于该试件整体的综合渗透率值，是一个单一方向上的参数。它不能或难以获得同一个试件在不同方向和不同测定尺度上的渗透率参数，更难以获得该试件渗透率场或空间分布的信息。

然而，大量的研究已经表明，岩石渗透率的空间分布是高度非均质的，并随测定样品尺度的变化而变化，为此在某一固定的样品尺度上获得的渗透率参数很难在其它尺度上应用。同时，岩石渗透率还是高度各向异性的，在某一单一方向上获得的渗透率也难以推广到其它方向上应用。因此，迫切需要一种能够针对低渗透岩芯特点可在多尺度和多个方向上测定岩石渗透率参数的突破性技术，以解决本领域科学研究和工程应用中的瓶颈问题。正是在此技术背景下，本发明提出了基于高密度的电性参数监测条件下的水力学实验获得不同尺度和方向上的岩芯渗透率参数的技术，为研究岩石渗透率的空间分布、尺度变化和各向异性特征提供了重要手段，可在各种生产工程的岩芯渗透率测定中广泛应用。相关技术，国内外均还未见报道。

本发明中所指的岩芯不仅是通过钻孔采芯所获得的岩石样品，而是更广泛意义上的圆柱体形状或多棱柱形状的多孔介质，其直径在厘米级至米级范围，高度也在厘米级至米级范围。因此，除岩石样品之外，还包括土壤、混凝土、木材等凡是能让流体渗透的多孔介质。本发明中所指的流体包括液体和气体，在液体中包括无极溶液和有机液体，如水、盐水、油类等凡是能够渗透到多孔介质中并引起介质电性差异的液体，气体中包括空气、氮气、二氧化碳等凡是能够渗透到多孔介质中的气体。本发明中的流体包括不同温度条件下的流体。本发明中所指的电性参数包括充电率值或极化率值、自然电位值以及不同测定频率下的电阻值、阻抗值和电容值。本发明中所指的渗透（或水力）实验是让流体进入岩芯或从岩芯中释出的实验。

技术实现要素：

通过在岩芯周围的高密度电极布置实现对实验过程中岩芯电性参数高密度监测的思想与方法（图1）。可用于监测的电性参数包括电阻值、电容值、充电率值和自然电位值。通过南京九州勘探技术有限公司研制的高密度电阻值测定系统ert21实现对任意电极间的电阻、电容、充电率值和自然电位值的动态扫描测定，获得岩芯中气体或液体的渗透过程，绘制出岩芯不同位置处电性参数的时间变化曲线（穿透曲线）（图2和3）。针对围压条件下的高密度电性参数测定，我们发明了“高密度电性参数监测室内专用岩芯夹持器”（正在申请专利）专门用于在围岩条件下的岩芯电性参数测定，实现了在不同实验条件下对岩芯内部电性变化过程的监测。实验结果表明，无论是液体还是气体其在岩芯内部的穿透过程都能够通过岩芯电性参数的变化体现出来。

通过对电性参数时间变化曲线的分析，提取出岩芯中渗透气体或液体到达时间的方法（图2和3）。根据电性参数穿透曲线，采用数学分析方法，依据对拐点、极值等的求取，获得各监测点处的流体到达时间。根据相邻监测点间的流体到达时间差异，获得流体在岩芯中的渗透速度。并由此获得流体在岩芯中到达时间和渗透速度的三维空间分布结构。由于电极分布在岩芯周围不同的高度和方向上，因此所得的到达时间和速度代表着不同的测定尺度和不同的方向，从而实现了在同一岩芯上不同尺度和方向上的参数测定。

针对流体在岩芯中渗透过程的基于四面体微元剖分的三维数值模拟方法（图4和5）。采用所研发的软件系统reticu3d（已申请软件著作权登记）和tetramesh（已申请软件著作权登记）对所测定的岩芯进行三维四面体剖分，结合实验条件建立起模拟所需的边界和初始条件，采用改进的tough2软件系统对电极监测处的流体饱和度变化过程进行模拟监测，获得模拟的流体到达时间。

通过对模拟流体到达时间与全部实测流体到达时间的拟合，获得岩芯综合渗透率参数的方法（图6和7）。为了消除在测定和模拟中流体到达绝对时间误差的影响，本发明采用流体到达各电极之间的时间差进行拟合，提高了拟合结果的可靠性。由于模拟结果还会受到样品初始含水量以及水分特征曲线中有关参数的影响，我们还通过对各参数的敏感度分析，找出各参数的选择范围，最终确定最佳渗透率参数。

通过模拟流体到达时间与不同方向和测定尺度上获得的流体到达时间的分别拟合，获得同一个岩芯在不同方向上和不同测定尺度上的渗透率参数的方法。设置在岩芯不同高度和方向上的电极监测了流体在该方向和尺度上的流动过程，对其流体到达时间的拟合即可获得该方向和尺度下的渗透率参数，从而在同一个岩芯上方便的实现了多尺度和多方向的渗透率测定。

通过实测流体到达的时间场，建立起对前述获得的渗透率参数的修正系数场，实现非均匀渗透率场条件下的模拟到达时间与实测到达时间拟合的思想与方法（图8、9、10、11和12）。在用样的条件下，到达时间早的区域，流体渗透速度快，其渗透率值也较大，反之则渗透率值小。通过对实测到达时间的归一化，建立起全域的渗透率修正系数场。在tough2系统中进行非均匀渗透率场条件下的模拟，拟合模拟获得的到达时间差与实测到达时间差，最终获得实验岩芯渗透率参数在空间上的分布。

附图说明

图1为电极在岩芯周围高密度布置和电性参数测定方法示意图，其中左图是针对自然电位的监测方案，右图是针对电阻、充电率以及电容参数的监测方案。所布置的电极数量可以根据要求增减，监测也参数可以选择和组合。由于每两两电极之间就可以获得一个电性参数，一次扫描测定就可以获得几十甚至上百个的电性参数值。扫描监测动态进行，每几分钟就可以获得一次扫描数据，再结合传统的试件两端的压力监测，从而实现对岩芯渗透过程的多参数全程监测。

图2和3分别是流体渗透过程中e号和b号岩芯电阻值随时间的变化过程以及流体到达监测点时间的确定，其中电阻值的大小作了归一化处理。对于花岗岩之类高电阻率的岩体，在实验初期会因为岩芯导电性较差而出现数据的不稳定现象，这是正常的。流体到达时间我们在曲线的最速下降段求取，但也可以用曲线的拐点位置来定义，因为我们在后面的参数求取中，是以监测点之间的时间差异来计算的，到达时间的系统增减并不影响渗透率参数的最终计算结果。

图4是针对e号（左图）和b号（右图）两个岩芯的四面体剖分方案（仅显示剖分方案的表面）。对试件采用四面体剖分有利于准确刻画试件的不规则边界并在不同精度下实现三维条件下的模拟计算。为此，我们专门编制了tetramesh系统用于生成tough2所接收的剖分文件，实现了对所测岩芯的三维非饱和流计算。

图5是e号（左图）和b号（右图）岩芯上测定位置处模拟的微元饱和度随时间变化曲线以及流体到达时间的提取。流体的到达以饱和度逐渐增加的方式出现，我们在曲线的最速上升段上确定模拟的流体到达时间。

图6和7分别是针对e号和b号岩芯监测微元处模拟的流体到达时间与实测时间的比较。小图(a)是各测定位置处实际到达时间与模拟到达时间的比较，小图(b)是对二者相关关系的分析。小图(c)是监测点间的实测时间差与模拟到达的时间差的比较，小图(d)是两种时间差的相关关系分析，理想情况是数据点聚集在1:1线附近。

图8分别是针对e号（左图）和b号（右图）岩芯的渗透率修正系数在各监测点上的分布。颜色越深气泡越大表示修正系数越大。渗透率修正系数范围限定在0.3至1.8之间。

图9是采用非均匀的渗透率场进行模拟后，监测点处的流体饱和度随时间的变化曲线以及所确定的流体到达时间，左图为e号岩芯，右图为b号岩芯。与图5相比，可以看出，各监测点的流体到达时间差异明显。

图10和11分别是e号和b号岩芯采用非均质渗透率场模拟后，模拟到达时间与实测到达时间的比较。可以清晰的看出，采用非均质渗透率场模拟后，模拟与实测的到达时间能够很好的吻合。

图12分别是针对e号和b号岩芯所获得的渗透率场的空间分布。上层指上面一圈测定电极所在的截面，下层则是下面一圈电极所在的截面。

具体实施方式

（一）进行高密度电性参数监测条件下的水力实验

如图1所示，根据测定要求，在岩芯周围高密度布置电极，在岩芯水力实验的同时，通过所布置的电极，进行电极间自然电位、电阻、电容、充电率等电性参数的扫描监测。监测持续到入渗流体到达甚至超过下层电极设置的位置为止。

（二）提取流体到达时间

根据所获得的电性参数穿透曲线，提取入渗水体到达监测位置处的时间（图2和3），通过监测位置之间流体到达的时间差求取流体入渗的速度，确定监测位置的座标，建立起流体到达时间或速度的文件。

（三）针对岩芯生成四面体有限元剖分网格

对岩芯进行三维四面体剖分，是一项具有许多技巧并很难保证一次成功的工作。对圆柱形的岩芯进行三维四面体剖分的方法较多，本发明采用我们研发的镜像剖分方法。其基本思想是将圆柱体水平分割为多个薄层，对其中的一个薄层进行四面体剖分，然后再通过镜像映射获得与该薄层镜像对称的薄层的剖分方案，最后再将各个部分组合起来，获得整个圆柱体的最终剖分方案（图4）。如果通过镜像映射即可获得整个圆柱体的剖分，则无需后续组合即实现剖分，这在圆柱体的层数是2的整数次方时适用。如果剖分圆柱体的层数不是2的整数次方，则通过与奇数层剖分结果的组合实现对整个圆柱体的四面体有限元剖分。由于镜像法将复杂的目标体分解为单层的简单目标体，减少了剖分节点数量，因而很容易实现对该单层的剖分，从而保证了对整个圆柱体的完全剖分。此外，为了使剖分的网格具有一定的密度同时又保证剖分的成功，本发明先采用较粗的网格进行剖分，然后再通过所研发的tetradetail命令对网格进行加密。由于该命令可以反复使用，从而可以获得针对圆柱体的高密度四面体剖分网格。

（四）建立起模拟所需的边界和初始条件

本发明采用tough2系统对岩芯中的渗透过程进行三维模拟，当然也可以采用其它软件系统对岩芯中的渗透过程进行三维模拟。对于tough2系统，采用具有大体积（e+50）和小距离（e-09）的微元来定义给定水头边界，也可以采用静默微元来定义给定水头边界，这都可以通过我们所研发的软件命令tetramesh、blockmesh或iblockmesh来实现。对于定通量和变通量边界则可以采用所研发的generation命令来实现。其中的定水头边界值（单位：pa），根据实验中施加的压力再加上大气压力计算而得，例如实验中施加的压力为0.55mpa，则定水头边界的压力值为0.55mpa+1.013e+5pa=6.513e+5pa。水头按每厘米98.0665pa计算。对于初始条件，可以通过模拟时的参数输入文件来定义。

（五）确定模拟监测点位置提取出监测处的微元序列

根据前述流体到达时间文件中的监测点座标，在岩芯四面体剖分网格中提取出相应的微元名称，构成模拟监测微元序列，以便模拟相应微元中的饱和度随时间的变化过程。

（六）建立模拟用的输入参数文件，将监测用的微元序列放入到输入参数文件中，修改渗透率和相关模拟参数，进行三维的流体入渗过程模拟

用于tough2模拟模块eos9输入的参数文件内容和说明见后（其中有汉字的行是说明行，实际执行时需将该行删除），可在linux下通过命令vim对其进行修改。初始条件采用外部文件incon给定，此时给定的为定水头边界方式。如果需要采用通量边界，则给定定义通量边界的文件gener。

rocks----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8物名行数颗粒密度孔隙度渗透率kx渗透率ky渗透率kz导热系数比热容grnte22600.e00.027.e-178.e-177.e-172.51920.k类型lambda值饱和度sr饱和度ss饱和度sgl7.257.101.0.10p类型lambda值饱和度sr1/参照p0最大pmax饱和度ss7.257.050.9869e-51.e71.参照条件参照压力参照温度水面参照高程强制用srrefco1.013e515.0.2start----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8允许incon中定义的微元顺序是随机的----*----1-mop:123456789*123456789*1234----*----5----*----6----*----7----*----8param----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8步长控制解法总步间隔123456789012345-7890-2342500501010000002000040003开始时间终止时间额外步长重力加速度08.64e6-1.9.81步长1步长21.32e-41.32e4收敛标准1.e-5s初始值0.05foft----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8所要模拟监测的微元名序列a1d52a1c23a1976a1a25a1c95a1e9a1x88a1w59a1u13a1u61a1x32a1y45times----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8在特定时间输出（每行个数总数）11时间表1.e-9indom----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8按区域进行初始条件定义incon----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8按微元进行初始条件定义，可用单独文件定义此时无需本数据块微元名孔隙度a1a120.2000000e-01主变量1初始值主变量2初始值1.1012598500000e+050.0000000000000e+00gener----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8源汇项定义，可用单独文件定义此时无需本数据块nover----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8禁止版本信息等输出endcy----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8文件结束并开始模拟。

（七）比较监测微元处流体的到达时间与测定获得的到达时间，调整渗透率等参数，反复模拟直至两种方法所获得的流体到达时间吻合，由此确定岩芯的（综合）渗透率参数

对监测微元处流体到达时间的提取通过对输出结果文件foft的分析实现，即绘制饱和度随时间的变化曲线，提取各监测微元处入渗水体的到达时间（图5）。比较模拟到达时间和实测到达时间（图6和7），调整计算参数。重复上述计算，直至所得到达时间差和实测到达时间差的相关性达到可接受的程度，即认为模拟成功，模拟中所采用的参数即为测定样品的参数。在实例中，由此获得e号岩芯的参数为：孔隙度n=0.07，渗透率k=2.0e-17m²。b号岩芯的参数为：孔隙度n=0.10，渗透率k=6.0e-17m²。通过与文献报道结果及其它测试方法获得结果的比较可知，此处所得的渗透率参数是比较合理的。

(八)对同一岩芯中不同方向上的渗透率参数的计算

显然，布置在岩芯不同方向上的电极所监测到的流体是来自于该方向上的流体，对该方向上流体到达时间的拟合所得到的渗透率即是该方向上的岩芯渗透率。因此实验测定时，应在同一高度上尽可能在不同的位置和方向上布置较多的电极，实现对多个方向的监测。同时还需要在模拟中采用严格对称的剖分网格，以免由于剖分网格的不对称性，造成计算结果在方向上的人为差异。在此基础上，针对岩芯不同方向上的渗透率可以通过两种拟合方式获得：一是分别针对不同方向上的实测到达时间进行拟合，获得不同方向上的渗透率。二是针对所有方向上的实测到达时间进行拟合，通过调整x、y和z三个方向上的渗透率参数，来获得最佳拟合结果。因此，本发明提供了一种在同一个岩芯上方便的测定其不同方向上的渗透率参数的技术，避免了传统方法采用不同的岩芯获得不同方向上的渗透率时所可能引起的误差。

(九)对同一岩芯中不同尺度下的渗透率参数的计算

由于不同高度上的电极所监测到的流体是流经不同路径长度后到达的流体，因此对不同高度上所获得的流体到达时间分别进行拟合，即可获得不同岩芯尺度条件下的渗透率参数。例如，在实例中，e号岩样的上层电极距顶面距离为l=16mm，所获得的渗透率为kx=2.0e-17m²,ky=2.0e-17m²,kz=2.0e-17m²，下层电极距离顶面的距离为l=36mm，所获得的渗透率为kx=0.6e-17m²,ky=8.0e-17m²,kz=2.0e-17m²。对于b号岩样，其上层电极l=25mm，所得渗透率为kx=2.0e-17m²,ky=4.0e-17m²,kz=4.0e-17m²，而其下层电极l=52mm，相应的渗透率值为kx=4.0e-17m²,ky=11e-17m²,kz=6.0e-17m²。由于上层电极监测的流体所流过的路径显然要小于下层电极监测水体流过的路径，因此根据上层监测数据获得的是小尺度条件下的渗透率参数，而下层获得的则是稍大些尺度条件下的渗透率参数。通过在不同高度上多层电极的布置和监测，即可实现对多尺度渗透率参数的测定。根据e号和b号岩样的结果，渗透率具有随着测定尺度增加而增大的趋势。由此可见，本发明提供了一种在同一个岩芯上即可方便的实现多尺度渗透率测定的技术，避免了采用多个不同尺度的岩芯获得不同尺度条件下的渗透率所引起的误差。

（十）对同一岩芯中非均质渗透率场的计算

通过水力实验条件下的高密度电性参数监测，获得岩芯渗透率的空间分布结构，是本发明的一个重要特色，在测试体的渗透率参数空间分布结构识别中具有非常重要的意义。本发明的基本思路是，通过对所测定的流体到达时间场的分析，建立起渗透率修正系数的空间分布结构，再根据三维空间上修正后的渗透率场模拟流体的入渗过程，使其与实测的流体到达时间完美吻合。

渗透率修正系数场的建立通过归一化实测的流体到达时间，再通过转换和插值来获得（图8）。在此基础上，通过所研发的软件命令tetramesh建立起含渗透率修正系数的tough2模拟用输入文件mesh，执行含渗透率修正系数的三维数值模拟，拟合监测微元处的流体到达时间与实测到达时间。调整修正系数范围，反复模拟，直至获得最佳拟合结果。例如，图9给出了非均质渗透率场条件下e号和b号岩芯样品模拟的监测微元饱和度的时间变化过程，显然渗透率的非均质性使得入渗水体到达同样高度的时间有了明显的差别，尤其是对于下层电极，因为水体入渗途径更长，其到达时间的差异也就更大，这和实测结果是完全一致的。图10和11是非均质渗透率条件下e号和b号岩芯的模拟到达时间与实测时间的比较，从中可以看出监测微元处流体到达时间的吻合程度得到了很大的提高，说明介质渗透率的非均质性确实是影响模拟结果的重要因素。同时，也证明根据入渗水体的到达时间场来获得非均质渗透率场的方法是可行和可靠的。

根据渗透率修正系数结合前述所得的综合渗透率参数即可获得岩芯不同位置的渗透率参数，实现对岩芯非均质渗透率场的获取。显然，监测电极的数量越多，该渗透率场与实际吻合的程度就越高。图12给出了e号（左图）和b号（右图）岩芯在其上层和下层电极平面上获得的渗透率空间分布，岩芯渗透率空间分布的非均匀性显而易见。