CN109284555A - 基于多孔介质几何形状估计渗透率的网格随机游走方法 - Google Patents

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

Info

Publication number
CN109284555A
CN109284555A CN201811135485.XA CN201811135485A CN109284555A CN 109284555 A CN109284555 A CN 109284555A CN 201811135485 A CN201811135485 A CN 201811135485A CN 109284555 A CN109284555 A CN 109284555A
Authority
CN
China
Prior art keywords
random walk
permeability
grid
porous media
penetration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201811135485.XA
Other languages
English (en)
Other versions
CN109284555B (zh
Inventor
南统超
吴吉春
蒋建国
曾献奎
王远坤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University
Original Assignee
Nanjing University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University filed Critical Nanjing University
Priority to CN201811135485.XA priority Critical patent/CN109284555B/zh
Publication of CN109284555A publication Critical patent/CN109284555A/zh
Application granted granted Critical
Publication of CN109284555B publication Critical patent/CN109284555B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/08Probabilistic or stochastic CAD

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于多孔介质几何形状估计渗透率的网格随机游走方法,该方法首先将孔隙空间从数值化的多孔介质中分离并离散化;假定孔隙与固体骨架的交界面为吸收边界,介质外壁为反弹边界;选定介质外壁上一点,用行走器在孔隙离散网格上重复模拟若干次随机游走实现,行走器一旦接触吸收边界立即被吸收,导致行走终止;对出发点集内每个出发点都重复若干次随机游走;行走器被吸收的位置相对于外壁的深度为穿透深度,对所有的穿透深度做平均,得到估计的平均穿透深度;根据平均穿透深度估算多孔介质渗透率。与现有技术相比,本发明方法精度更高,而且计算简便、效率高,对孔隙形状无特殊要求,适用性更强。

Description

基于多孔介质几何形状估计渗透率的网格随机游走方法
技术领域
本发明属于水力学技术领域,具体涉及一种基于多孔介质几何形状估计渗透率的网格随机游走方法。
背景技术
多孔介质渗透率是地下水、岩石物理、油气开采等领域常用的基础性参数,在地下水开采、地下地表水交互、农业灌溉管理、放射性废物处置等问题中都有重要的科学与工程意义。作为岩土体的宏观参数,渗透率与孔隙微结构密切相关。尽管借助传统方法在实验室中通过渗透实验可以测量多孔介质渗透系数和渗透率,这类方法对岩土体样品采取与制备、测量仪器硬件、实验操作等要求都比较高,而且耗时长。随着新型探测手段如核磁共振、X射线显微成像、扫描电镜成像等技术的应用,多孔介质的微观几何结构数据收集正变得日渐简便、普遍。根据孔隙几何特征直接估算介质渗透率成为了可能。Hwang等人于2000年在“On the rapid estimation of permeability for porous media using Brownianmotion paths”一文中提出孔隙内部随机游走模拟的平均穿透深度D可以体现渗透率k的大小,并取k=D2。Simonov与Mascagni在2004年的“Random Walk Algorithms forestimating Effective Properties of Digitized Porous Media”一文中考虑了孔隙度n的影响,认为k=nD2,并且使用了“球面随机游走”和“立方体随机游走”两种方式模拟孔隙中的随机游走。然而,“球面随机游走”算法在判断当前点到最近边界的距离时比较困难、费时,而且很难准确判断抵达边界的时间、位置。Nan与Wu于2018年在“Random Walk Pathsolution to groundwater flow dynamics in highly heterogeneous aquifers”一文中认为“立方体随机游走”算法只能应用于内外边界都是多面体的孔隙,对形状卷曲、粗糙、不规则的孔隙以及尖灭的狭缝等无法直接应用,必须经过几何形状的多面体简化,由此造成额外误差。Sabelfeld在2013年于“A stochastic spectral projection method forsolving PDEs in domains composed by overlapping discs,spheres,and half-spaces”一文中提出一种“谱映射”方法,对于固体骨架可看作部分重叠的圆盘或球体的多孔介质,可以用近似解表示而避免数值计算的困难,但对粗糙、不规则性强的多孔介质并不适用。此外,经检验,以往的平均穿透深度计算渗透率的公式k=D2和k=nD2带有一定偏差,需要进行校正。
发明内容
针对于上述现有技术的不足,本发明的目的在于提供一种基于多孔介质几何形状估计渗透率的网格随机游走方法,该方法将网格随机游走与平均穿透深度估计法有机结合,以解决现有技术中对内部粗糙、不规则的一般性多孔介质无法通过现有“球面随机游走”和“立方体随机游走”等方法准确估计渗透率的问题。
为达到上述目的,本发明采用的技术方案如下:
本发明的一种基于多孔介质几何形状估计渗透率的网格随机游走方法,包括步骤如下:
(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=nCD2估计渗透率,其中,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“随机游走循环”),如此就完成了一个网格随机游走实现。Monte Carlo模拟要求对同一个出发点重复模拟N次随机游走实现(图2中“Monte Carlo循环”),记录每次的穿透深度di,j(参照图1,No.8;图2步骤(f))。从出发点集内每个可能的出发点出发都进行N次随机游走模拟(参照图2步骤(g)和“出发点循环”),记录所有的穿透深度di,j,并计算所有di,j值的平均值,得到平均穿透深度D(图2步骤(h))。利用改进的渗透率估计公式k=nCD2估算介质渗透率(图2步骤(i)),其中C≈1.125,n为孔隙度。通过对不同形状孔隙的渗透率的计算,发现本发明的方法与解析解吻合的很好;和现有随机游走估计渗透率的方法相比,本发明获得的精度更高;对于界面形状弯曲、带有尖灭的孔隙空间,本发明也可有效应用并获得足够精确的计算结果。
下面结合具体实施例对本发明做进一步的解释,其中涉及一些简写符号,以下为注解:
D:平均穿透深度;
k:渗透率;
N:Monte Carlo模拟次数;
WOG:本发明中的网格随机游走法;
WOS:球体随机游走方法;
WOC:立方体随机游走法。
实施例1:Hele-Shaw模型
研究区为介于两平行平板之间的空间(参照图3),平板间距为l,若l足够小,板间水流为Stokes流。此模型的孔隙度n=1.0,渗透率有解析解:k=l2/12=0.0833…l2
采用WOG,将板间空间用l/100的步长离散并进行网格随机游走,取Monte Carlo模拟次数N=105,获得的平均穿透深度D≈0.276l,估算的渗透率kG=0.0859l2。WOS估计的渗透率为0.0709l2。WOC估计的渗透率为0.0763l2
实施例2:横截面为边长l的正方形的直孔隙
研究对象为柱状直孔隙,横截面是边长为l的正方形(参照图4),l足够小。此模型的孔隙度n=1.0,渗透率有解析解:k≈0.0351l2
采用WOG方法估计渗透率,将孔隙空间用l/100的步长离散并进行网格随机游走,取Monte Carlo模拟次数N=105,获得的平均穿透深度D≈0.175l,估算的渗透率kG=0.0345l2。WOS估计的渗透率为0.0307l2。WOC估计的渗透率为0.0289l2
实施例3:横截面为半径R的圆形的直孔隙
研究对象为柱状直孔隙,横截面是半径为R的圆形(参照图5),R足够小。此模型的孔隙度n=1.0,渗透率有解析解:k=0.125R2
采用WOG方法估计渗透率,将孔隙空间用R/100的步长离散并进行网格随机游走,取Monte Carlo模拟次数N=105,获得的平均穿透深度D≈0.3318R,估算的渗透率kG=0.124R2。WOS估计的渗透率为0.105R2。WOC不适用于曲面形状,无计算结果。
实施例4:等大的固体球间孔隙空间
研究对象为等大的固体球整齐排列后留下的孔隙空间(参照图6),固体球直径为l=1mm。此模型孔隙度为0.476,渗透率没有解析解。为了有参考结果可以比较,使用数值模拟软件COMSOL模拟该介质内的Stokes流,并根据流量计算得渗透率约为2.52×10-3l2。采用WOG方法估计渗透率,网格离散步长为l/100,模拟次数N=104,所有吸收点位置参见图7(坐标轴单位为l/100)。获得的平均穿透深度D≈0.0759l,估算的渗透率为kG=3.09×10-3l2。而WOS估计的平均穿透深度0.0960l,估计的渗透率为4.39×10-3l2。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更高。
本发明具体应用途径很多,以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进,这些改进也应视为本发明的保护范围。

Claims (6)

1.一种基于多孔介质几何形状估计渗透率的网格随机游走方法,其特征在于,包括步骤如下:
(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估计渗透率。
2.根据权利要求1所述的基于多孔介质几何形状估计渗透率的网格随机游走方法,其特征在于,所述步骤(1)中的将孔隙空间用网格进行剖分离散,具体包括:采用有限元法或有限差分法的离散方式,基于均匀网格或非均匀网格对孔隙空间进行剖分离散。
3.根据权利要求1所述的基于多孔介质几何形状估计渗透率的网格随机游走方法,其特征在于,所述步骤(2)中的吸收边界会将到达该处的行走器吸收,使该次随机游走终止;反弹边界将到达该处的行走器反弹回孔隙空间。
4.根据权利要求1所述的基于多孔介质几何形状估计渗透率的网格随机游走方法,其特征在于,所述步骤(3)具体包括:当前网格点到相邻网格点的转移概率通过相关点之间的相对距离计算得到。
5.根据权利要求1所述的基于多孔介质几何形状估计渗透率的网格随机游走方法,其特征在于,所述步骤(4)中在介质外壁上出发点集的选定按照均匀分布随机抽取,或者包含每个可能的网格点。
6.根据权利要求1所述的基于多孔介质几何形状估计渗透率的网格随机游走方法,其特征在于,所述步骤(8)具体包括:通过平均穿透深度D及公式k=nCD2估计渗透率,其中,k为渗透率,n为孔隙度,C≈1.125。
CN201811135485.XA 2018-09-28 2018-09-28 基于多孔介质几何形状估计渗透率的网格随机游走方法 Active CN109284555B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811135485.XA CN109284555B (zh) 2018-09-28 2018-09-28 基于多孔介质几何形状估计渗透率的网格随机游走方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811135485.XA CN109284555B (zh) 2018-09-28 2018-09-28 基于多孔介质几何形状估计渗透率的网格随机游走方法

Publications (2)

Publication Number Publication Date
CN109284555A true CN109284555A (zh) 2019-01-29
CN109284555B CN109284555B (zh) 2023-02-14

Family

ID=65181665

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811135485.XA Active CN109284555B (zh) 2018-09-28 2018-09-28 基于多孔介质几何形状估计渗透率的网格随机游走方法

Country Status (1)

Country Link
CN (1) CN109284555B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109918748A (zh) * 2019-02-25 2019-06-21 南京大学 基于随机游走的非均质含水层水流问题格林函数计算方法
CN111026028A (zh) * 2019-12-11 2020-04-17 上海维宏电子科技股份有限公司 数控系统中针对加工工件实现二维平面化网格划分处理的方法
CN111563927A (zh) * 2020-05-14 2020-08-21 西南石油大学 一种基于岩石微ct图像的孔隙迂曲度计算方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106574981A (zh) * 2015-08-17 2017-04-19 数岩科技(厦门)股份有限公司 针对多孔介质的核磁共振分析系统和方法
CN107808049A (zh) * 2017-10-26 2018-03-16 南京大学 基于多孔介质三维微观结构模型的dnapl迁移数值模拟方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106574981A (zh) * 2015-08-17 2017-04-19 数岩科技(厦门)股份有限公司 针对多孔介质的核磁共振分析系统和方法
US20170122858A1 (en) * 2015-08-17 2017-05-04 Irock Technologies Co., Ltd. Nmr anaylysis system and method for porous media
CN107808049A (zh) * 2017-10-26 2018-03-16 南京大学 基于多孔介质三维微观结构模型的dnapl迁移数值模拟方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
金毅等: "煤微观结构三维表征及其孔-渗时空演化模式数值分析", 《岩石力学与工程学报》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109918748A (zh) * 2019-02-25 2019-06-21 南京大学 基于随机游走的非均质含水层水流问题格林函数计算方法
CN111026028A (zh) * 2019-12-11 2020-04-17 上海维宏电子科技股份有限公司 数控系统中针对加工工件实现二维平面化网格划分处理的方法
CN111026028B (zh) * 2019-12-11 2022-12-06 上海维宏电子科技股份有限公司 针对加工工件实现二维平面化网格划分处理的方法
CN111563927A (zh) * 2020-05-14 2020-08-21 西南石油大学 一种基于岩石微ct图像的孔隙迂曲度计算方法

Also Published As

Publication number Publication date
CN109284555B (zh) 2023-02-14

Similar Documents

Publication Publication Date Title
Counillon et al. Flow-dependent assimilation of sea surface temperature in isopycnal coordinates with the Norwegian Climate Prediction Model
Yeh et al. An iterative stochastic inverse method: Conditional effective transmissivity and hydraulic head fields
Moulton et al. The black box multigrid numerical homogenization algorithm
CN109284555A (zh) 基于多孔介质几何形状估计渗透率的网格随机游走方法
Liu et al. A multilevel sampling method for detecting sources in a stratified ocean waveguide
Li et al. Strengthened linear sampling method with a reference ball
CN107766666B (zh) 一种基于分数阶差分法的三维时域电磁反常扩散模拟方法
Sargsyan et al. Uncertainty quantification given discontinuous model response and a limited number of model runs
Wołoszyn et al. Sensitivity analysis of efficiency thermal energy storage on selected rock mass and grout parameters using design of experiment method
CN109632604B (zh) 一种孔隙尺度到岩心尺度聚合物驱相对渗透率粗化方法
Zhang et al. Modeling free-surface seepage flow in complicated fractured rock mass using a coupled RPIM-FEM method
Zhao et al. Stochastic modeling of the permeability of randomly generated porous media via the Lattice Boltzmann method and probabilistic collocation method
Hanke One shot inverse scattering via rational approximation
CN109582922B (zh) 一种水合物饱和度的现场快速判别方法及处理终端
Porta et al. A space–time adaptation scheme for unsteady shallow water problems
Arbogast et al. Multiscale mortar mixed methods for heterogeneous elliptic problems
Shahrokhabadi et al. Random isogeometric analysis for modeling seepage in unsaturated soils
Ghosh et al. Effect of initial compaction state on near-saturated hydraulic conductivity
Mahdavi et al. A localized differential quadrature model for moving boundary shallow water flows
CN114428043A (zh) 多孔介质孔径分布表征方法及电子设备
Wang et al. Characteristics of spatial distribution of soil water-air-heat parameters in typical oasis croplands at middle reaches of Heihe River
CN114137623B (zh) 核磁测井仪观测模式确定方法、存储介质以及电子设备
Bevilacqua et al. Enhancing the uncertainty quantification of pyroclastic density current dynamics in the Campi Flegrei caldera.
Liu et al. Applying COSISIM model to study the permeability of porous media
CN109918748B (zh) 基于随机游走的非均质含水层水流问题评估方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant