CN111797562B - 裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备 - Google Patents

裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备 Download PDF

Info

Publication number
CN111797562B
CN111797562B CN202010697717.1A CN202010697717A CN111797562B CN 111797562 B CN111797562 B CN 111797562B CN 202010697717 A CN202010697717 A CN 202010697717A CN 111797562 B CN111797562 B CN 111797562B
Authority
CN
China
Prior art keywords
crack
fracture
representing
determining
characteristic unit
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.)
Active
Application number
CN202010697717.1A
Other languages
English (en)
Other versions
CN111797562A (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.)
Tsinghua University
Original Assignee
Tsinghua 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 Tsinghua University filed Critical Tsinghua University
Priority to CN202010697717.1A priority Critical patent/CN111797562B/zh
Publication of CN111797562A publication Critical patent/CN111797562A/zh
Application granted granted Critical
Publication of CN111797562B publication Critical patent/CN111797562B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备。该方法包括：S100：构建三维离散裂隙网络模型；S200：确定该模型的截面上的裂隙开度分布，并确定截面上每个特征单元的裂隙初始开度；S300：确定每个特征单元的水压力，并水压力确定水流速；S400：确定每个特征单元新的裂隙开度；S500：比较每个特征单元的新的裂隙开度和与其对应的裂隙初始开度，当至少一个特征单元的新的裂隙开度不满足预设条件时，将新的裂隙开度作为该特征单元的新的裂隙初始开度，返回执行S300，当每个特征单元的新的裂隙开度均满足预设条件时，执行S600；S600：确定三维裂隙岩体的溶质浓度分布。能够为深部地下岩石工程成功建设和安全运行提供重要的理论基础。

Description

裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备
技术领域
本发明涉及裂隙岩体渗流与溶质运移的数值模拟领域，尤其涉及一种裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备。
背景技术
近年来，随着地下资源的开发和利用，高放废物地质处置、CO2地质封存、页岩气开采以及增强型地热系统等重要工程均涉及三维裂隙岩体中的渗流和物质运移问题。在裂隙岩体中，原岩的孔隙率通常比较小(<1％)，故其渗透性通常可忽略，而裂隙岩体的渗流及溶质传输通道则主要由分布在其中的裂隙网络所提供。受人为或自然扰动的影响，地应力重分布会导致裂隙开度显著变化，进而对裂隙岩体中的地下水渗流和溶质传输过程产生影响。所以，全面深入的理解力学过程影响下的溶质(如放射性核素)传输过程是深部地下岩石工程成功建设和安全运行的重要理论基础。
相较于孔隙岩体，裂隙岩体具有高度的非均质性，这导致对裂隙岩体中力学 -渗流耦合过程的数值模拟面临巨大困难。目前，用于描述裂隙岩体中力学-渗流耦合作用的模型主要包括等效连续体模型和离散裂隙网络模型。相对于等效连续体模型，离散裂隙网络模型可以更真实地刻画单裂隙的局部行为以及复杂裂隙网络的宏观行为，并能很好地刻画由于裂隙随机分布和开度不均引起的渗流各向异性和非均质性，并因此得到了广泛的应用。
粒子追踪方法，即随机行走模拟方法，将溶质视为大量粒子的集合，这些粒子的分布即可视为溶质的分布，该方法通过模拟大量粒子的运动来模拟溶质运移过程。粒子追踪方法是模拟裂隙岩体中流体与溶质运移的常用方法。其适用范围比较广泛，且非常适合模拟裂隙岩体中的溶质运移过程。
裂隙岩体中的力学-渗流耦合作用是通过应力对裂隙的法向压闭、张开和剪胀作用，显著改变裂隙的开度和水力传导特性，进而影响渗流场。反过来，渗流场通过渗透压力影响裂隙的应力场和变形。因此，在裂隙岩体溶质运移问题中，需要考虑应力和渗流相互作用下的裂隙中的平流和基质扩散过程。现有技术中的方案只针对孔隙介质或单条裂隙等简单情况，且并未考虑力学-渗流的耦合作用，缺乏对复杂裂隙网络尤其是三维裂隙网络中溶质运移过程的研究。
发明内容
本发明的主要目的是提供一种裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备，以解决三维裂隙岩体中考虑力学-渗流耦合作用下预测三维裂隙岩体溶质浓度分布的问题。
第一方面，本申请的实施方式提供一种三维裂隙岩体溶质浓度分布的预测方法，包括以下步骤：S100：获取三维裂隙岩体的裂隙网络几何统计参数，根据裂隙网络几何统计参数构建三维离散裂隙网络模型；S200：截取三维离散裂隙网络模型的截面来表征三维裂隙岩体的裂隙面，利用网格将所述截面划分为若干个特征单元，基于傅立叶积分法确定三维离散裂隙网络模型的截面上的裂隙开度分布，并根据所述截面上的裂隙开度分布确定所述截面上每个特征单元的裂隙初始开度；S300：根据立方定律构建三维裂隙岩体的渗流模型，根据每个特征单元的裂隙初始开度，利用三维裂隙岩体的渗流模型确定每个特征单元的水压力，并根据每个特征单元的水压力确定每个特征单元的水流速；S400：对于每个特征单元，根据裂隙初始开度和水压力，基于刚块弹簧法确定新的裂隙开度；S500：比较每个特征单元的新的裂隙开度和与其对应的裂隙初始开度，当至少一个特征单元的新的裂隙开度和与其对应的裂隙初始开度的差值的绝对值大于或等于预设阈值时，将新的裂隙开度作为该特征单元的新的裂隙初始开度，返回执行S300，当每个特征单元的新的裂隙开度和与其对应的裂隙初始开度的差值的绝对值均小于预设阈值时，执行S600；S600：基于粒子追踪法，确定基质岩石中的溶质浓度，并根据基质岩石中的溶质浓度以及每个特征单元的水流速和新的裂隙开度，确定三维裂隙岩体的溶质浓度分布。
在一个实施例中，基于傅立叶积分法确定三维离散裂隙网络模型的截面上的裂隙开度分布，包括：确定三维离散裂隙网络模型的截面上的裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数；对所述指数协方差函数进行采样，确定采样协方差序列；对采样协方差序列进行傅立叶变换，得到谱密度函数；根据谱密度函数，确定振幅谱；确定相位谱；根据振幅谱和相位谱，确定复傅立叶系数；对复傅立叶系数进行傅立叶逆变换，得到三维离散裂隙网络模型的截面上的裂隙开度分布。
在一个实施例中，裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数的表达式为：
Figure GDA0002750390750000031
其中，C(k)表示裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数，σ2表示方差，λ表示相关长度，|k|表示所述截面上两特征单元的指定顶点之间的距离；
采样协方差序列的表达式为：
{C(k)；k＝1,…,N}
其中，N表示采样数，k表示第k个特征单元；
谱密度函数的表达式为：
Figure GDA0002750390750000032
其中，S(j)表示离散谱密度函数，j表示频率；
振幅谱的表达式为：
Figure GDA0002750390750000033
其中，A(j)表示振幅谱；
相位谱的表达式为：
Figure GDA0002750390750000034
其中，
Figure GDA0002750390750000035
表示相位谱，Uj表示均匀分布于[0,1]的随机变量；
复傅立叶系数的表达式为：
Figure GDA0002750390750000036
其中，
Figure GDA0002750390750000037
表示复傅立叶系数，i表示复数单位，R(j)表示实部，为偶数，I(j)表示虚部，为奇数；
三维裂隙岩体的裂隙开度分布的表达式为：
Figure GDA0002750390750000041
其中，z(k)表示第k个特征单元的裂隙初始开度。
在一个实施例中，确定每个特征单元的水压力，包括：
利用下式确定每个特征单元的水压力：
Figure GDA0002750390750000042
其中，pm和pn分别表示特征单元m和特征单元n的水压力，Qmn表示特征单元m和特征单元n之间的流量，
Figure GDA0002750390750000043
Rmn表示特征单元m 和特征单元n之间的流阻，υ表示水的动力粘度，Δx和Δy分别表示特征单元的长度和宽度，bm和bn分别表示特征单元m和特征单元n的裂隙初始开度。
在一个实施例中，根据每个特征单元的水压力确定该特征单元的水流速，包括：
利用下式确定每个特征单元的水流速：
Figure GDA0002750390750000044
其中，vm表示特征单元m的水流速，υ表示水的动力粘度，
Figure GDA0002750390750000045
表示特征单元m的压力梯度，bm表示特征单元m的裂隙初始开度。
在一个实施例中，对于每个特征单元，根据裂隙初始开度和水压力，基于刚块弹簧法确定新的裂隙开度，包括：基于简化Barton-Bandis模型，根据每个特征单元的水压力确定每个特征单元中裂隙的法向位移；基于Coulomb滑移准则，确定每个特征单元中裂隙的剪胀量；根据每个特征单元的裂隙初始开度、裂隙的法向位移和裂隙的剪胀量，确定该特征单元的新的裂隙开度。
在一个实施例中，利用下式确定每个特征单元中裂隙的法向位移：
un′m＝(σn′m-pn′m)/kn′m
其中，σn′m表示特征单元m中的裂隙法向应力，un′m表示特征单元m中裂隙的法向位移，pn′m表示特征单元m的水压力，kn′m表示裂隙的法向刚度。
在一个实施例中，基于Coulomb滑移准则，确定每个特征单元中裂隙的剪胀量，包括：根据每个特征单元所承受的切向应力，确定每个特征单元中裂隙的切向位移；根据每个特征单元中裂隙的切向位移，确定每个特征单元中裂隙的剪胀量。
在一个实施例中，利用下式确定每个特征单元中裂隙的切向位移：
σsm＝usmks
其中，σsm表示特征单元m中裂隙所承受的切向应力，usm表示特征单元m 中裂隙的切向位移，ks表示裂隙的切向刚度。
在一个实施例中，利用下式确定每个特征单元中裂隙的剪胀量：
Figure GDA0002750390750000051
其中，udm表示特征单元m中裂隙的剪胀量，usm表示特征单元m中裂隙的切向位移，up表示裂隙的峰值切向位移，ucs表示裂隙的峰后切向位移。
在一个实施例中，根据每个特征单元的裂隙初始开度、裂隙的法向位移和裂隙的剪胀量，确定该特征单元的新的裂隙开度，包括：
利用下式确定该特征单元的新的裂隙开度：
bcm＝bm+un′m+udm
其中，bcm表示力学-渗流耦合作用下特征单元m的新的裂隙开度，bm表示特征单元m的裂隙初始开度，un′m为特征单元m的裂隙的法向位移，udm表示特征单元m的裂隙的剪胀量。
在一个实施例中，基于粒子追踪法，根据每个特征单元的水流速和新的裂隙开度，确定三维裂隙岩体的溶质浓度分布，包括：
基于溶质粒子在基质岩石中的扩散传输机理，确定基质岩石中的溶质浓度，包括：利用下式确定基质岩石中的溶质浓度：
Figure GDA0002750390750000052
其中，Cr表示基质岩石中的溶质浓度，Dr表示基质岩石的扩散系数，
Figure GDA0002750390750000053
表示梯度算子；
基于溶质粒子在裂隙中的平流传输机理，根据基质岩石中的溶质浓度以及每个特征单元的水流速和新的裂隙开度，确定裂隙中的溶质浓度，包括：利用下式确定裂隙中的溶质浓度：
Figure GDA0002750390750000061
其中，Cfm表示特征单元m中裂隙的溶质浓度，t表示时间，ufm表示特征单元m的水流速，
Figure GDA0002750390750000062
表示所述裂隙面上的梯度算子，Df表示裂隙中的扩散系数， Cr表示基质岩石中的溶质浓度，αs表示质量传递系数，bcm表示力学-渗流耦合作用下特征单元m的新的裂隙开度。
第二方面，本申请的实施方式提供一种存储介质，存储有计算机程序，所述计算机程序被处理器执行时，实现如上文所述的三维裂隙岩体溶质浓度分布的预测方法的步骤。
第三方面，本申请的实施方式提供一种计算机设备，包括处理器和存储有程序代码的存储介质，所述程序代码被所述处理器执行时，实现如上文所述的三维裂隙岩体溶质浓度分布的预测方法的步骤。
根据现场测得的裂隙分布参数，包括：裂隙迹长、间距、密度、开度和走向等信息，建立的三维离散裂隙网络模型，同时考虑了裂隙开度的非均质分布，相较于单裂隙或二维离散裂隙网络模型更接近岩体和裂隙的真实分布；基于粒子追踪方法对三维裂隙岩体中的溶质运移进行模拟，采用粒子模拟溶质(放射性核素)，通过计算粒子在力学-渗流耦合作用下的运动轨迹(随时间的位置变化)，在统计意义上反映溶质在裂隙网络中的浓度和迁移速度等参数，并考虑溶质在裂隙网络中的平流和基质扩散等传输机理，可以更精确地刻画溶质(放射性核素)在离散裂隙网络中的运移。为深部地下岩石工程成功建设和安全运行提供重要的理论基础。
附图说明
构成本申请的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定，在附图中：
图1为根据本申请一示例性实施方式的三维裂隙岩体溶质浓度分布的预测方法的流程图；
图2为根据本申请一具体实施例的三维离散裂隙网络模型的示意图；
图3为根据本申请一具体实施例的水流过三维离散裂隙网络模型的截面的示意图；
图4为根据本申请一具体实施例的三维离散裂隙网络模型的截面上的特征单元的分布示意图；
图5A为基于刚块弹簧法确定裂隙的法向位移的理论曲线图；
图5B为基于刚块弹簧法确定裂隙的剪胀量的理论曲线图。
具体实施方式
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。
实施例一
如图1所示，本申请的实施方式提供一种三维裂隙岩体溶质浓度分布的预测方法，包括以下步骤：
S100：获取三维裂隙岩体的裂隙网络几何统计参数，根据裂隙网络几何统计参数构建三维离散裂隙网络模型。
其中，裂隙网络几何统计参数可以通过现场测量得到，其可以包括：裂隙迹长、间距、密度和走向等信息。图2为根据本申请一具体实施例的三维离散裂隙网络模型的示意图。如图2所示，对于所建立的三维离散裂隙网络模型，将相邻岩块的界面定义为裂隙。
S200：截取三维离散裂隙网络模型的截面来表征三维裂隙岩体的裂隙面，利用网格将所述截面划分为若干个特征单元，基于傅立叶积分法确定三维离散裂隙网络模型的截面上的裂隙开度分布，并根据所述截面上的裂隙开度分布确定所述截面上每个特征单元的裂隙初始开度。
其中，用来将三维离散裂隙网络模型的截面划分为若干个特征单元的网格可以为矩形网格或三角形网格等。
基于傅立叶积分法确定三维离散裂隙网络模型的截面上的裂隙开度分布，可以包括以下步骤：
第一步，确定三维离散裂隙网络模型的截面上的裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数。具体的，裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数的表达式可以为：
Figure GDA0002750390750000081
其中，C(k)表示裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数，σ2表示方差，λ表示相关长度，|k|表示所述截面上两特征单元的指定顶点之间的距离。
第二步，对所述指数协方差函数进行采样，确定采样协方差序列。具体的，采样协方差序列的表达式可以为：
{C(k)；k＝1,…,N}
其中，N表示采样数，k表示第k个特征单元。
第三步，对采样协方差序列进行傅立叶变换，得到谱密度函数。离散谱密度函数的表达式可以为：
Figure GDA0002750390750000082
其中，S(j)表示离散谱密度函数，j表示频率。
第四步，根据谱密度函数，确定振幅谱。振幅谱的表达式可以为：
Figure GDA0002750390750000083
其中，A(j)表示离散振幅谱。
第五步，确定相位谱。相位谱的表达式可以为：
Figure GDA0002750390750000084
其中，
Figure GDA0002750390750000085
表示相位谱，Uj表示均匀分布于[0,1]的随机变量；
第六步，根据振幅谱和相位谱，确定复傅立叶系数。实函数的傅立叶变换是厄米函数，复傅立叶系数的表达式可以为：
Figure GDA0002750390750000091
其中，
Figure GDA0002750390750000092
表示复傅立叶系数，i表示复数单位，R(j)表示实部，为偶数，I(j)表示虚部，为奇数。
第七步，对复傅立叶系数
Figure GDA0002750390750000093
进行傅立叶逆变换，得到第二步中的采样协方差序列的空间实现，即裂隙开度的非均质分布，也即三维离散裂隙网络模型的截面上的裂隙开度分布。三维裂隙岩体的裂隙开度分布的表达式可以为：
Figure GDA0002750390750000094
其中，z(k)表示第k个特征单元的裂隙初始开度。
S300：根据立方定律构建三维裂隙岩体的渗流模型，根据每个特征单元的裂隙初始开度，利用三维裂隙岩体的渗流模型确定每个特征单元的水压力，并根据每个特征单元的水压力确定每个特征单元的水流速。
具体的，图3为根据本申请一具体实施例的水流过三维离散裂隙网络模型的截面的示意图。如图3所示，非均质的裂隙开度分布包含了若干个特征单元，将立方定律应用到每一个特征单元，两个相邻特征单元，其裂隙开度分别为bm和 bn，当从特征单元m到特征单元n的流量为Qmn时，两个单元的压力差可写为：
Figure GDA0002750390750000095
图4为根据本申请一具体实施例的三维离散裂隙网络模型的截面上的特征单元的分布示意图。如图4所示，根据质量守恒方程，通过每个特征单元m的净流量为0，即：可以利用下式确定每个特征单元的水压力(下式中的n求和不是指对所有的特征单元求和，而是指对与特征单元m邻接的特征单元进行水流量求和，对于特征单元m来说，流出量与流入量之差为0)：
Figure GDA0002750390750000096
其中，pm和pn分别表示特征单元m和特征单元n的水压力，Qmn表示特征单元m和特征单元n之间的流量，
Figure GDA0002750390750000097
Rmn表示特征单元m 和特征单元n之间的流阻，υ表示水的动力粘度，Δx和Δy分别表示特征单元的长度和宽度，bm和bn分别表示特征单元m和特征单元n的裂隙初始开度。
具体的，为每一个特征单元应用上述质量平衡方程，形成一个稀疏线性方程组，即：
Figure GDA0002750390750000101
该方程组包含了N2个未知数和N2个方程(定解方程组)，通过迭代法求解可得到裂隙网络中的水压力分布。
进一步的，根据每个特征单元的水压力确定该特征单元的水流速，可以包括：利用下式确定每个特征单元的水流速：
Figure GDA0002750390750000102
其中，vm表示特征单元m的水流速，υ表示水的动力粘度，
Figure GDA0002750390750000103
表示特征单元m的压力梯度，bm表示特征单元m的裂隙初始开度。
S400：对于每个特征单元，根据裂隙初始开度和水压力，基于刚块弹簧法确定新的裂隙开度，可以包括：
第一步，基于简化Barton-Bandis模型，根据每个特征单元的水压力确定每个特征单元中裂隙的法向位移。图5A为基于刚块弹簧法确定裂隙的法向位移的理论曲线图。
当法向应力σn′增加(加载)时，法向刚度kn′呈分段式非线性增大(如图5A 所示，法向刚度在图中表现为斜率)，法向应力σn′1、σn′2和σn′3分别对应法向位移un′1、un′2和un′3，法向应力当法向应力减小(卸载)时，法向刚度kn′与应力加载曲线重合。因此，在某一法向应力σn′m的作用下，考虑裂隙水压力pm，特征单元m中裂隙的法向位移可以表示为：
un′m＝(σn′m-pn′m)/kn′m
其中，σn′n表示特征单元m中的裂隙法向应力，un′m表示特征单元m中裂隙的法向位移，pn′m表示特征单元m的水压力，kn′m表示裂隙的法向刚度。
第二步，基于Coulomb滑移准则，确定每个特征单元中裂隙的剪胀量，可以包括：根据每个特征单元所承受的切向应力，确定每个特征单元中裂隙的切向位移；根据每个特征单元中裂隙的切向位移，确定每个特征单元中裂隙的剪胀量。
图5B为基于刚块弹簧法确定裂隙的剪胀量的理论曲线图。
当切向应力低于剪切强度时，剪切刚度ks保持恒定(如图5B所示)，当切向应力高于剪切强度时，切向应力保持恒定，此时剪胀量ud开始线性增加，剪胀角为φd，而当切向位移us大于ucs时，剪胀停止。对于特征单元m，可以利用下式确定每个特征单元中裂隙的切向位移：
σsm＝usmks
其中，σsm表示特征单元m中裂隙所承受的切向应力，usm表示特征单元m 中裂隙的切向位移，ks表示裂隙的切向刚度。
可以利用下式确定每个特征单元中裂隙的剪胀量：
Figure GDA0002750390750000111
其中，udm表示特征单元m中裂隙的剪胀量，usm表示特征单元m中裂隙的切向位移，up表示裂隙的峰值切向位移，ucs表示裂隙的峰后切向位移。
第三步，根据每个特征单元的裂隙初始开度、裂隙的法向位移和裂隙的剪胀量，确定该特征单元的新的裂隙开度：可以包括：利用下式确定该特征单元的新的裂隙开度：
bcm＝bm+un′m+udm
其中，bcm表示力学-渗流耦合作用下特征单元m的新的裂隙开度，bm表示特征单元m的裂隙初始开度，un′m为特征单元m的裂隙的法向位移，udm表示特征单元m的裂隙的剪胀量。
当计算出整个裂隙面上所有特征单元的新的裂隙开度后，再重新用于渗流分析，直到变形和渗流同时满足收敛条件，即小于一定的容差，也即执行S500。
S500：比较每个特征单元的新的裂隙开度和与其对应的裂隙初始开度，
当至少一个特征单元的新的裂隙开度和与其对应的裂隙初始开度的差值的绝对值大于或等于预设阈值时，将新的裂隙开度作为该特征单元的新的裂隙初始开度，返回执行S300，
当每个特征单元的新的裂隙开度和与其对应的裂隙初始开度的差值的绝对值均小于预设阈值时，执行S600。
S600：基于粒子追踪法，确定基质岩石中的溶质浓度，并根据基质岩石中的溶质浓度以及每个特征单元的水流速和新的裂隙开度，确定三维裂隙岩体的溶质浓度分布。
以浓度大小表征粒子的运移过程，裂隙岩体的溶质传输过程涉及到裂隙内部物质的对流扩散以及与基质岩石的物质交换，到达交叉点的溶质粒子沿某一裂隙流出的概率等于该裂隙的流量与流出该交叉点的总流量之比。
具体的，基于溶质粒子在基质岩石中的扩散传输机理，确定基质岩石中的溶质浓度，可以包括：利用下式的控制方程确定基质岩石中的溶质浓度：
Figure GDA0002750390750000121
其中，Cr表示基质岩石中的溶质浓度，Dr表示基质岩石的扩散系数，
Figure GDA0002750390750000122
表示梯度算子。
基于溶质粒子在裂隙中的平流传输机理，根据基质岩石中的溶质浓度以及每个特征单元的水流速和新的裂隙开度，确定裂隙中的溶质浓度，可以包括：利用下式的控制方程确定裂隙中的溶质浓度：
Figure GDA0002750390750000123
其中，Cfm表示特征单元m中裂隙的溶质浓度，t表示时间，ufm表示特征单元m的水流速，
Figure GDA0002750390750000131
表示所述裂隙面上的梯度算子，Df表示裂隙中的扩散系数， Cr表示基质岩石中的溶质浓度，αs表示质量传递系数，bcm表示力学-渗流耦合作用下特征单元m的新的裂隙开度。
通过迭代法求解上述两个控制方程即可得到三维裂隙岩体的溶质浓度分布。通过统计流出边界上的溶质粒子浓度，可以得到溶质粒子浓度随时间的变化规律。
根据现场测得的裂隙分布参数，包括：裂隙迹长、间距、密度、开度和走向等信息，建立的三维离散裂隙网络模型，同时考虑了裂隙开度的非均质分布，相较于单裂隙或二维离散裂隙网络模型更接近岩体和裂隙的真实分布；基于粒子追踪方法对三维裂隙岩体中的溶质运移进行模拟，采用粒子模拟溶质(放射性核素)，通过计算粒子在力学-渗流耦合作用下的运动轨迹(随时间的位置变化)，在统计意义上反映溶质在裂隙网络中的浓度和迁移速度等参数，并考虑溶质在裂隙网络中的平流和基质扩散等传输机理，可以更精确地刻画溶质(放射性核素)在离散裂隙网络中的运移。本发明能够为深部地下岩石工程成功建设和安全运行提供重要的理论基础。
实施例二
本申请的实施方式提供一种存储介质，存储有计算机程序，所述计算机程序被处理器执行时，实现如上文所述的三维裂隙岩体溶质浓度分布的预测方法的步骤。
实施例三
本申请的实施方式提供一种计算机设备，包括处理器和存储有程序代码的存储介质，所述程序代码被所述处理器执行时，实现如上文所述的三维裂隙岩体溶质浓度分布的预测方法的步骤。
需要说明的是，本文可提供包含特定值的参数的示范，但这些参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应值。
还需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。
除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
应当理解的是，本说明书中的示例性实施方式可以由多种不同的形式来实施，并且不应当被解释为只限于这里所阐述的实施方式。提供这些实施方式是为了使得本申请的公开彻底且完整，并且将这些示例性实施方式的构思充分传达给本领域普通技术人员，而不应当理解为对本发明的限制。

Claims (7)

1.一种三维裂隙岩体溶质浓度分布的预测方法，其特征在于，包括以下步骤：
S100：获取三维裂隙岩体的裂隙网络几何统计参数，根据裂隙网络几何统计参数构建三维离散裂隙网络模型；
S200：截取三维离散裂隙网络模型的截面来表征三维裂隙岩体的裂隙面，利用网格将所述截面划分为若干个特征单元，基于傅立叶积分法确定三维离散裂隙网络模型的截面上的裂隙开度分布，并根据所述截面上的裂隙开度分布确定所述截面上每个特征单元的裂隙初始开度；
S300：根据立方定律构建三维裂隙岩体渗流模型，根据每个特征单元的裂隙初始开度，利用三维裂隙岩体渗流模型确定每个特征单元的水压力，并根据每个特征单元的水压力确定每个特征单元的水流速；
S400：对于每个特征单元，根据裂隙初始开度和水压力，基于刚块弹簧法确定新的裂隙开度，包括：
基于简化Barton-Bandis模型，根据每个特征单元的水压力确定每个特征单元中裂隙的法向位移，包括：利用下式确定每个特征单元中裂隙的法向位移：
un′m＝(σn′m-pn′m)/kn′m
其中，σn′m表示特征单元m中的裂隙法向应力，un′m表示特征单元m中裂隙的法向位移，pn′m表示特征单元m的水压力，kn′m表示裂隙的法向刚度；
基于Coulomb滑移准则，确定每个特征单元中裂隙的剪胀量，包括：根据每个特征单元所承受的切向应力，确定每个特征单元中裂隙的切向位移；根据每个特征单元中裂隙的切向位移，确定每个特征单元中裂隙的剪胀量；其中，利用下式确定每个特征单元中裂隙的切向位移：
σsm＝usmks
其中，σsm表示特征单元m中裂隙所承受的切向应力，usm表示特征单元m中裂隙的切向位移，ks表示裂隙的切向刚度；
利用下式确定每个特征单元中裂隙的剪胀量：
Figure FDA0003986293010000021
其中，udm表示特征单元m中裂隙的剪胀量，usm表示特征单元m中裂隙的切向位移，up表示裂隙的峰值切向位移，ucs表示裂隙的峰后切向位移；
根据每个特征单元的裂隙初始开度、裂隙的法向位移和裂隙的剪胀量，确定该特征单元的新的裂隙开度，包括：利用下式确定该特征单元的新的裂隙开度：
bcm＝bm+un′m+udm
其中，bcm表示力学-渗流耦合作用下特征单元m的新的裂隙开度，bm表示特征单元m的裂隙初始开度，un′m为特征单元m的裂隙的法向位移，udm表示特征单元m的裂隙的剪胀量；
S500：比较每个特征单元的新的裂隙开度和与其对应的裂隙初始开度，
当至少一个特征单元的新的裂隙开度和与其对应的裂隙初始开度的差值的绝对值大于或等于预设阈值时，将新的裂隙开度作为该特征单元的新的裂隙初始开度，返回执行S300，
当每个特征单元的新的裂隙开度和与其对应的裂隙初始开度的差值的绝对值均小于预设阈值时，执行S600；
S600：基于粒子追踪法，确定基质岩石中的溶质浓度，并根据基质岩石中的溶质浓度以及每个特征单元的水流速和新的裂隙开度，确定三维裂隙岩体的溶质浓度分布，包括：
基于溶质粒子在基质岩石中的扩散传输机理，确定基质岩石中的溶质浓度，包括：利用下式确定基质岩石中的溶质浓度：
Figure FDA0003986293010000022
其中，Cr表示基质岩石中的溶质浓度，Dr表示基质岩石的扩散系数，
Figure FDA0003986293010000023
表示梯度算子；
基于溶质粒子在裂隙中的平流传输机理，根据基质岩石中的溶质浓度以及每个特征单元的水流速和新的裂隙开度，确定裂隙中的溶质浓度，包括：利用下式确定裂隙中的溶质浓度：
Figure FDA0003986293010000031
其中，Cfm表示特征单元m中裂隙的溶质浓度，t表示时间，ufm表示特征单元m的水流速，
Figure FDA0003986293010000032
表示所述裂隙面上的梯度算子，Df表示裂隙中的扩散系数，Cr表示基质岩石中的溶质浓度，αs表示质量传递系数，bcm表示力学-渗流耦合作用下特征单元m的新的裂隙开度。
2.根据权利要求1所述的三维裂隙岩体溶质浓度分布的预测方法，其特征在于，基于傅立叶积分法确定三维离散裂隙网络模型的截面上的裂隙开度分布，包括：
确定三维离散裂隙网络模型的截面上的裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数；
对所述指数协方差函数进行采样，确定采样协方差序列；
对采样协方差序列进行傅立叶变换，得到谱密度函数；
根据谱密度函数，确定振幅谱；
确定相位谱；
根据振幅谱和相位谱，确定复傅立叶系数；
对复傅立叶系数进行傅立叶逆变换，得到三维离散裂隙网络模型的截面上的裂隙开度分布。
3.根据权利要求2所述的三维裂隙岩体溶质浓度分布的预测方法，其特征在于，裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数的表达式为：
Figure FDA0003986293010000033
其中，C(k)表示裂隙初始开度与特征单元的指定顶点的空间坐标之间的指数协方差函数，σ2表示方差，λ表示相关长度，|k|表示所述截面上两特征单元的指定顶点之间的距离；
采样协方差序列的表达式为：
{C(k)；k＝1，...，N}
其中，N表示采样数，k表示第k个特征单元；
谱密度函数的表达式为：
Figure FDA0003986293010000041
其中，S(j)表示谱密度函数，j表示频率；
振幅谱的表达式为：
Figure FDA0003986293010000042
其中，A(j)表示振幅谱；
相位谱的表达式为：
Figure FDA0003986293010000043
其中，
Figure FDA0003986293010000044
表示相位谱，Uj表示均匀分布于[0，1]的随机变量；
复傅立叶系数的表达式为：
Figure FDA0003986293010000045
其中，
Figure FDA0003986293010000046
表示复傅立叶系数，i表示复数单位，R(j)表示实部，为偶数，I(j)表示虚部，为奇数；
三维裂隙岩体的裂隙开度分布的表达式为：
Figure FDA0003986293010000047
其中，z(k)表示第k个特征单元的裂隙初始开度。
4.根据权利要求1所述的三维裂隙岩体溶质浓度分布的预测方法，其特征在于，确定每个特征单元的水压力，包括：
利用下式确定每个特征单元的水压力：
Figure FDA0003986293010000048
其中，pm和pn分别表示特征单元m和特征单元n的水压力，Qmn表示特征单元m和特征单元n之间的流量，
Figure FDA0003986293010000051
Rmn表示特征单元m和特征单元n之间的流阻，υ表示水的动力粘度，Δx和Δy分别表示特征单元的长度和宽度，bm和bn分别表示特征单元m和特征单元n的裂隙初始开度。
5.根据权利要求1所述的三维裂隙岩体溶质浓度分布的预测方法，其特征在于，根据每个特征单元的水压力确定该特征单元的水流速，包括：
利用下式确定每个特征单元的水流速：
Figure FDA0003986293010000052
其中，vm表示特征单元m的水流速，υ表示水的动力粘度，
Figure FDA0003986293010000053
表示特征单元m的压力梯度，bm表示特征单元m的裂隙初始开度。
6.一种存储介质，存储有计算机程序，其特征在于，所述计算机程序被处理器执行时，实现如权利要求1-5中任一项所述的三维裂隙岩体溶质浓度分布的预测方法的步骤。
7.一种计算机设备，包括处理器和存储有程序代码的存储介质，所述程序代码被所述处理器执行时，实现如权利要求1-5中任一项所述的三维裂隙岩体溶质浓度分布的预测方法的步骤。
CN202010697717.1A 2020-07-20 2020-07-20 裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备 Active CN111797562B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010697717.1A CN111797562B (zh) 2020-07-20 2020-07-20 裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010697717.1A CN111797562B (zh) 2020-07-20 2020-07-20 裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备

Publications (2)

Publication Number Publication Date
CN111797562A CN111797562A (zh) 2020-10-20
CN111797562B true CN111797562B (zh) 2023-03-24

Family

ID=72807962

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010697717.1A Active CN111797562B (zh) 2020-07-20 2020-07-20 裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备

Country Status (1)

Country Link
CN (1) CN111797562B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106815460A (zh) * 2016-10-11 2017-06-09 中国辐射防护研究院 一种离散裂隙网络评价放射性核素在岩石裂隙中迁移的方法
CN106886682A (zh) * 2017-01-04 2017-06-23 中国环境科学研究院 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法
CN109887083A (zh) * 2019-01-29 2019-06-14 中国石油集团测井有限公司西南分公司 一种裂缝-孔隙双重介质耦合渗透率模型的建立方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR3041026B1 (fr) * 2015-09-15 2017-10-20 Ifp Energies Now Procede pour caracteriser le reseau de fractures d'un gisement fracture et procede pour l'exploiter

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106815460A (zh) * 2016-10-11 2017-06-09 中国辐射防护研究院 一种离散裂隙网络评价放射性核素在岩石裂隙中迁移的方法
CN106886682A (zh) * 2017-01-04 2017-06-23 中国环境科学研究院 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法
CN109887083A (zh) * 2019-01-29 2019-06-14 中国石油集团测井有限公司西南分公司 一种裂缝-孔隙双重介质耦合渗透率模型的建立方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Effects of Fracture Surface Roughness on Macroscopic Fluid Flow and Solute Transport in Fracture Networks;Zhihong Zhao等;《Rock Mesh Rock Eng》;20131109;全文 *
三维双重孔隙-裂隙介质热-水-应力-迁移耦合模型及其有限元分析;张玉军等;《中国科学:技术科学》;20101215(第12期);全文 *
压力溶解和自由面溶解/沉淀对双重孔隙岩体中热-水-应力耦合影响的有限元分析;张玉军等;《岩石力学与工程学报》;20111215(第12期);全文 *
岩体裂隙中渗流场有限元随机模拟分析;李守巨等;《岩土力学》;20090710(第07期);全文 *

Also Published As

Publication number Publication date
CN111797562A (zh) 2020-10-20

Similar Documents

Publication Publication Date Title
Hunt et al. Flow, transport, and reaction in porous media: Percolation scaling, critical‐path analysis, and effective medium approximation
Grenier et al. Groundwater flow and heat transport for systems undergoing freeze-thaw: Intercomparison of numerical simulators for 2D test cases
Lei et al. The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks
Flühler et al. Lateral solute mixing processes—a key for understanding field-scale transport of water and solutes
Haldorsen Simulator parameter assignment and the problem of scale in reservoir engineering
US6230101B1 (en) Simulation method and apparatus
Raymond et al. Numerical analysis of thermal response tests with a groundwater flow and heat transfer model
US8078437B2 (en) Upscaling reservoir models by reusing flow solutions from geologic models
Rodriguez et al. High-resolution numerical simulation of flow through a highly sinuous river reach
Spanoudaki et al. Development and verification of a 3-D integrated surface water–groundwater model
US20150338550A1 (en) Method and system for characterising subsurface reservoirs
CN106837315B (zh) 裂缝性碳酸盐岩基质与裂缝耦合作用表征方法
Diamantopoulos et al. Inverse modeling of dynamic nonequilibrium in water flow with an effective approach
Lee et al. Fracture-based modeling of complex flow and CO2 migration in three-dimensional fractured rocks
Shein Physically based mathematical models in soil science: history, current state, problems, and outlook (analytical review)
CN103477248A (zh) 用于估算碳氢化合物生产区域中至少一个参数值的计算方法,以在该区域规划和实施作业
Bagtzoglou et al. Probabilistic simulation for reliable solute source identification in heterogeneous porous media
CN111797562B (zh) 裂隙岩体溶质浓度分布预测方法、存储介质和计算机设备
Suzuki et al. Characterization of 3D printed fracture networks
Vassilevski et al. CFD technology for 3D simulation of large-scale hydrodynamic events and disasters
CN114611434A (zh) 一种体现渣土改良泡沫剂渗流特性的模型构建方法及系统
Qiu et al. Numerical analysis of 1‐D compression wave propagation in saturated poroelastic media
Sun et al. Flow simulation in 3D fractured porous medium using a generalized pipe-based cell-centered finite volume model with local grid refinement
Patel et al. Direct numerical simulations of water flooding process through digitized porous rocks
Esmaeilpour Scale-Dependent petrophysical properties in porous media: A pore-network study

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