基于多孔介质几何形状估计渗透率的网格随机游走方法与流程

本发明属于水力学技术领域，具体涉及一种基于多孔介质几何形状估计渗透率的网格随机游走方法。

背景技术：

多孔介质渗透率是地下水、岩石物理、油气开采等领域常用的基础性参数，在地下水开采、地下地表水交互、农业灌溉管理、放射性废物处置等问题中都有重要的科学与工程意义。作为岩土体的宏观参数，渗透率与孔隙微结构密切相关。尽管借助传统方法在实验室中通过渗透实验可以测量多孔介质渗透系数和渗透率，这类方法对岩土体样品采取与制备、测量仪器硬件、实验操作等要求都比较高，而且耗时长。随着新型探测手段如核磁共振、x射线显微成像、扫描电镜成像等技术的应用，多孔介质的微观几何结构数据收集正变得日渐简便、普遍。根据孔隙几何特征直接估算介质渗透率成为了可能。hwang等人于2000年在“ontherapidestimationofpermeabilityforporousmediausingbrownianmotionpaths”一文中提出孔隙内部随机游走模拟的平均穿透深度d可以体现渗透率k的大小，并取k＝d²。simonov与mascagni在2004年的“randomwalkalgorithmsforestimatingeffectivepropertiesofdigitizedporousmedia”一文中考虑了孔隙度n的影响，认为k＝nd²，并且使用了“球面随机游走”和“立方体随机游走”两种方式模拟孔隙中的随机游走。然而，“球面随机游走”算法在判断当前点到最近边界的距离时比较困难、费时，而且很难准确判断抵达边界的时间、位置。nan与wu于2018年在“randomwalkpathsolutiontogroundwaterflowdynamicsinhighlyheterogeneousaquifers”一文中认为“立方体随机游走”算法只能应用于内外边界都是多面体的孔隙，对形状卷曲、粗糙、不规则的孔隙以及尖灭的狭缝等无法直接应用，必须经过几何形状的多面体简化，由此造成额外误差。sabelfeld在2013年于“astochasticspectralprojectionmethodforsolvingpdesindomainscomposedbyoverlappingdiscs,spheres,andhalf-spaces”一文中提出一种“谱映射”方法，对于固体骨架可看作部分重叠的圆盘或球体的多孔介质，可以用近似解表示而避免数值计算的困难，但对粗糙、不规则性强的多孔介质并不适用。此外，经检验，以往的平均穿透深度计算渗透率的公式k＝d²和k＝nd²带有一定偏差，需要进行校正。

技术实现要素：

针对于上述现有技术的不足，本发明的目的在于提供一种基于多孔介质几何形状估计渗透率的网格随机游走方法，该方法将网格随机游走与平均穿透深度估计法有机结合，以解决现有技术中对内部粗糙、不规则的一般性多孔介质无法通过现有“球面随机游走”和“立方体随机游走”等方法准确估计渗透率的问题。

为达到上述目的，本发明采用的技术方案如下：

本发明的一种基于多孔介质几何形状估计渗透率的网格随机游走方法，包括步骤如下：

(1)根据所要研究的多孔介质几何形状，将孔隙空间用网格进行剖分离散；

(2)孔隙-固体界面作为吸收边界，介质外壁作为反弹边界；

(3)根据网格剖分，计算每个网格点到相邻网格点的转移概率；

(4)在介质外壁上选取一个出发点集；

(5)以出发点集内的一出发点为起点，用网格随机游走法重复模拟n个随机游走，并记录n个穿透深度di，j；i为实现的序号，且i＝1，2，…，n；j为出发点序号；

(6)以出发点集内的下一个出发点为起点，同样重复模拟n个随机游走，记录穿透深度di，j+1，重复步骤(5)直至出发点集内所选定的出发点全部用完；

(7)对所有的穿透深度求平均值，得到平均穿透深度d；

(8)通过平均穿透深度d估计渗透率。

进一步地，所述步骤(1)中的将孔隙空间用网格进行剖分离散，具体包括：采用有限元法或有限差分法的离散方式，基于均匀网格或非均匀网格对孔隙空间进行剖分离散。均匀网格对称性好，离散误差更小；非均匀网格更适于局部网格加密；这两种剖分各有优势，在本发明中都能使用，可以根据具体问题灵活选择。

进一步地，所述步骤(2)中的吸收边界会将到达该处的行走器吸收，使该次随机游走终止；反弹边界将到达该处的行走器反弹回孔隙空间。

进一步地，所述步骤(3)具体包括：当前网格点到相邻网格点的转移概率通过相关点之间的相对距离计算得到；转移概率计算十分简单，而且每个网格点只需要在最初时计算一次并存储在计算机中。当行走器通过一个网格点，直接读取该网格点相关的转移概率。与其他随机游走方法相比，本发明方法仅涉及四则运算，计算量小，编程简便，算法执行效率高。

进一步地，所述步骤(4)中在介质外壁上出发点集的选定按照均匀分布随机抽取，或者包含每个可能的网格点，以保证所用的出发点集具有代表性。

进一步地，所述步骤(5)中，网格随机游走在步骤(1)的离散网格上进行而且仅从当前网格点移动到相邻网格点，虽然单次移动距离短，但单次移动计算量很小，而且可以直接判断是否到达吸收边界，而不必像其他方法那样每一步都花费额外时间去计算与边界的最小距离并判断是否到达边界。因此网格随机游走的计算效率比现有方法更高。

进一步地，所述步骤(6)中，出发点集内出发点的使用顺序不会影响最后的结果。即，某个网格点被作为第一个出发点或者是最后一个出发点，并不会对估计结果造成影响。

进一步地，所述步骤(8)具体包括：通过平均穿透深度d及公式k＝ncd²估计渗透率，其中，k为渗透率，n为孔隙度，c≈1.125，c为修正系数。现有的其他方法中默认取c＝1。本发明认为c＝1.125更合理，在所有存在解析解的模型中，c＝1.125的结果都比c＝1的结果更接近真实值，可以减小渗透率估计误差。

本发明的有益效果：

1.本发明可以将多孔介质渗透率估计问题转化为孔隙空间平均穿透深度的估计，不需渗流实验操作，避免实验成本；

2.与现有估计多孔介质渗透率的随机游走方法和谱映射法相比，本发明可以应用于界面粗糙、形状不规则的孔隙空间；

3.与现有估计多孔介质渗透率的随机游走方法相比，本发明不需要计算格林函数和边界距离，模拟更方便，操作更简单，效率也更高；

4.与现有估计多孔介质渗透率的随机游走方法相比，本发明估计的渗透率值更准确。

附图说明

图1为本发明中孔隙空间的网格随机游走模拟二维示意图。

图2为本发明方法于实施例中的执行流程图。

图3为由平行平板组成的hele-shaw模型中的网格随机游走实现示意图。

图4为横截面为边长l的正方形的直孔隙空间示意图。

图5为横截面为半径为r的圆形的直孔隙空间示意图。

图6为等大球形固体颗粒构成的多孔介质及其中的孔隙空间示意图。

图7为随机游走模拟出的二百一十万个吸收点位置示意图。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

本发明的基于多孔介质几何形状估计渗透率的网格随机游走方法，步骤如下：

图1中，no.1为介质固体骨架；no.2为孔隙-固体界面；no.3为介质外壁；no.4为孔隙内离散网格；no.5为当前选取的出发点；no.6为一个网格随机游走实现；no.7为随机游走实现的终结点，即行走器被吸收点；no.8为该次实现的穿透深度。

参照图1、图2所示，先构造一个孔隙空间上的离散网格(参照图1，no.4；图2步骤(a))，计算所有网格点到其相邻网格点的转移概率(图2步骤(b))。假设某网格点到六个方向(即x+、x-、y+，y-，z+和z-方向)上相邻点的距离分别为δx+、δx-、δy+、δy-、δz+、δz-，则到达x+方向邻点的概率可以按照距离反比加权法求得，即：px+＝1/δx+/(1/δx++1/δx-+1/δy++1/δy-+1/δz++1/δz-)；到其他相邻点的转移概率类似；若离散网格为均匀网格(即δx+＝δx-＝δy+＝δy-＝δz+＝δz-)，则各方向概率相等且均为1/6。

然后将介质外壁上网格点的集合作为出发点集，或者用均匀分布随机抽样从中获得一个子集作为出发点集，从出发点集选择一个出发点(图2步骤(c))，作为随机游走的起点，从该起点开始模拟随机游走(图2步骤(d))，根据转移概率随机移动行走器到相邻点，如此重复(图2步骤(e))直至行走器到达孔隙-固体界面被吸收(参照图1，no.7；图2“随机游走循环”)，如此就完成了一个网格随机游走实现。montecarlo模拟要求对同一个出发点重复模拟n次随机游走实现(图2中“montecarlo循环”)，记录每次的穿透深度di,j(参照图1，no.8；图2步骤(f))。从出发点集内每个可能的出发点出发都进行n次随机游走模拟(参照图2步骤(g)和“出发点循环”)，记录所有的穿透深度di,j，并计算所有di,j值的平均值，得到平均穿透深度d(图2步骤(h))。利用改进的渗透率估计公式k＝ncd²估算介质渗透率(图2步骤(i))，其中c≈1.125，n为孔隙度。通过对不同形状孔隙的渗透率的计算，发现本发明的方法与解析解吻合的很好；和现有随机游走估计渗透率的方法相比，本发明获得的精度更高；对于界面形状弯曲、带有尖灭的孔隙空间，本发明也可有效应用并获得足够精确的计算结果。

下面结合具体实施例对本发明做进一步的解释，其中涉及一些简写符号，以下为注解：

d：平均穿透深度；

k：渗透率；

n：montecarlo模拟次数；

wog：本发明中的网格随机游走法；

wos：球体随机游走方法；

woc：立方体随机游走法。

实施例1：hele-shaw模型

研究区为介于两平行平板之间的空间(参照图3)，平板间距为l，若l足够小，板间水流为stokes流。此模型的孔隙度n＝1.0，渗透率有解析解：k＝l²/12＝0.0833…l²。

采用wog，将板间空间用l/100的步长离散并进行网格随机游走，取montecarlo模拟次数n＝10⁵，获得的平均穿透深度d≈0.276l，估算的渗透率kg＝0.0859l²。wos估计的渗透率为0.0709l²。woc估计的渗透率为0.0763l²。

实施例2：横截面为边长l的正方形的直孔隙

研究对象为柱状直孔隙，横截面是边长为l的正方形(参照图4)，l足够小。此模型的孔隙度n＝1.0，渗透率有解析解：k≈0.0351l²。

采用wog方法估计渗透率，将孔隙空间用l/100的步长离散并进行网格随机游走，取montecarlo模拟次数n＝10⁵，获得的平均穿透深度d≈0.175l，估算的渗透率kg＝0.0345l²。wos估计的渗透率为0.0307l²。woc估计的渗透率为0.0289l²。

实施例3：横截面为半径r的圆形的直孔隙

研究对象为柱状直孔隙，横截面是半径为r的圆形(参照图5)，r足够小。此模型的孔隙度n＝1.0，渗透率有解析解：k＝0.125r²。

采用wog方法估计渗透率，将孔隙空间用r/100的步长离散并进行网格随机游走，取montecarlo模拟次数n＝10⁵，获得的平均穿透深度d≈0.3318r，估算的渗透率kg＝0.124r²。wos估计的渗透率为0.105r²。woc不适用于曲面形状，无计算结果。

实施例4：等大的固体球间孔隙空间

研究对象为等大的固体球整齐排列后留下的孔隙空间(参照图6)，固体球直径为l＝1mm。此模型孔隙度为0.476，渗透率没有解析解。为了有参考结果可以比较，使用数值模拟软件comsol模拟该介质内的stokes流，并根据流量计算得渗透率约为2.52×10^-3l²。采用wog方法估计渗透率，网格离散步长为l/100，模拟次数n＝10⁴，所有吸收点位置参见图7(坐标轴单位为l/100)。获得的平均穿透深度d≈0.0759l，估算的渗透率为kg＝3.09×10^-3l²。而wos估计的平均穿透深度0.0960l，估计的渗透率为4.39×10^-3l²。woc不适用于曲面形状，无计算结果。计算时间方面，在同一台pc机上，wog耗时约为2100秒；wos耗时约440000秒(边界厚度取l/100时)，时长约是前者的20倍。

以上四个实施例中的渗透率参考值及各方法的估计值、估计误差参照表1所示，表1如下：

表1

从上述表1可以看出，wog估计的渗透率精度要显著优于wos和woc。此外，woc仅对多面体空间适用；而wog的适用范围最广，与wos相当。wos虽然可以用于不规则孔隙空间，但计算边界距离要消耗大量计算时间，每次随机游走模拟时间约是wog的数倍以上。综合以上实施例的结果，可以得到如下结论：wog计算精度较现有方法有显著提高；算法比wos、woc简单；适用范围比woc更广；执行效率比wos更高。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进也应视为本发明的保护范围。