CN115661388A - 具有多尺度非均质性的三维裂隙网络岩体模型表征方法 - Google Patents

Abstract

本发明公开了一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，首先利用蒙特卡洛模拟方法生成随机三维离散裂隙网络岩体模型，对岩体模型进行三角网格和四面体网格划分，然后构建非均质性：赋予模型每条裂隙的平均开度值，即裂隙间宏观非均质性；赋予模型所有裂隙的平均开度值，即裂隙网络宏观非均质性；利用自仿射法和投影法赋予模型每条裂隙随机开度场，即裂隙内部开度微观尺度非均质性；最后对裂隙面三角网格单元和基质四面体网格单元进行数据转换，得到具有多尺度非均质性的三维裂隙网络岩体模型或网格模型。该方法能够更加真实地模拟复杂裂隙网络结构及周围岩石基质，为岩体裂隙渗流和污染物迁移的理论研究工作提供了准确的数值模型。

Description

具有多尺度非均质性的三维裂隙网络岩体模型表征方法

技术领域

本发明属于随机裂隙模型建立技术领域，尤其涉及一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法。

背景技术

岩体中大量发育的随机裂隙相互切割关系形成了空间连通的裂隙网络，这些连通性良好的裂隙网络构成了地下流体运动和污染物迁移的主要通道，对地下水渗流和污染物迁移以及核素迁移速率等特殊问题具有关键性控制作用。目前，表征复杂裂隙岩体的数学模型总体上分为等效连续介质模型、离散裂隙网络模型和裂隙网络岩体(或裂隙-孔隙双重介质)模型。与连续方法相比，离散裂隙网络方法因其准确考虑每条裂隙几何特征的优势而得到广泛应用，这种方法可以精确预测离散裂隙岩体内的流体流动和溶质迁移特征。

由于裂隙网络模型具有多尺度非均质性，准确表征裂隙岩体中溶质迁移过程仍具有一定的挑战性。该多尺度非均质性主要包括：(1)裂隙网络宏观尺度非均质性，如网络随机分布、裂隙密度和裂隙大小；(2)裂隙间开度分布宏观非均质性，即不同裂隙之间的裂隙开度具有变异性。一般地，不同裂隙的开度可以是相同值、随机值、与裂隙大小相关或半相关关系；(3)裂隙内部开度微观尺度非均质性，会造成每条裂隙内部的开度分布不同。近年来，研究人员对上述三种不同非均质尺度下的渗流(优势流通道、渗透率变化)和溶质迁移(“早到”和“拖尾”的非菲克迁移特征)开展了一些研究。例如，利用裂隙网络模型模型探讨裂隙随机分布和裂隙密度等网络尺度非均质性对不同溶质穿透曲线、平均击穿时间不确定性的影响。此外，一些学者的研究还发现，非均匀分布的裂隙开度会导致裂隙网络模型中形成极低流速和高流速区域，并形成显著的优先渗流路径，进而增强溶质的“早到”和“拖尾”特征。然而天然裂隙岩体由裂隙网络和多孔岩石组成，裂隙网络与周围岩石基质之间的相互作用对溶质迁移行为影响显著，特别是长时间尺度的溶质迁移过程。因此，考虑裂隙与岩石基质之间的物质交换过程对于准确预测裂隙岩体中的溶质迁移行为具有重要意义。近年来，随着计算机技术的发展与数值模拟技术的进步，对裂隙-基质系统中溶质迁移行为的显示模拟研究中，部分学者开始同时考虑三维裂隙网络分布和岩石基质的影响。

然而，由于三维离散裂隙网络几何特征较为复杂且具有随机分布的特点，导致裂隙网络与岩石基质的网格划分以及准确刻画和模拟两者之间的溶质交换仍有一定的难度，与之相关的数值模型或者模拟软件报道也较少。因此，构建同时考虑多尺度非均质性和岩石基质影响的三维裂隙网络岩体模型，有助于进一步深入研究复杂裂隙网络岩体中的流体流动和溶质迁移机制，对提高核废料处置、地下污染物治理等大型地下工程的安全稳定运行具有重要意义。

发明内容

针对上述问题，本发明的目的在于提供一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法。

为实现上述目的，本发明所采用的技术方案如下：

步骤1.利用蒙特卡洛模拟方法生成随机三维离散裂隙网络岩体模型；

步骤2.调用Delaunay三角网格开源工具箱DFNsMeshGenerator3D，对裂隙面和基质边界面进行三角网格划分，调用Delaunay四面体网格生成器TetGen对岩石基质进行四面体网格划分；

为了保证在进行裂隙面三角网格划分时，同时获得6个基质边界面的三角网格，在DFNsMeshGenerator3D输入文件中须增加6个裂隙面的顶点和中心坐标，将6个基质边界面进行三角网格划分；

为了保证在岩石基质四面体网格划分时，保持裂隙面及6个基质边界面的三角网格形态不变，TetGen的输入文件应包括裂隙面和基质边界面三角网格单元的节点列表和面列表；

步骤3.构建裂隙间宏观尺度非均质性或三维裂隙网络宏观尺度非均质性或裂隙内部开度微观尺度非均质性，并对随机三维离散裂隙网络岩体模型进行赋值；包括以下子步骤：

裂隙间宏观尺度非均质性：对模拟区域中的所有裂隙进行编号，根据裂隙开度与裂隙大小之间的幂律关系，计算得到每条裂隙的平均开度值b_correlated,i，对随机三维离散裂隙网络岩体模型每一条裂隙赋予开度值b_correlated,i，由于每条裂隙大小不同，由此得到的裂隙间开度差异即裂隙间宏观尺度非均质性；

裂隙网络宏观非均质性：根据随机三维离散裂隙网络岩体模型中每条裂隙的平均裂隙开度b_correlated,i，计算得到所有裂隙的平均开度值b_constant，对随机三维离散裂隙网络岩体模型每一条裂隙的平均裂隙开度值赋予b_constant，由此得到的裂隙网络差异即裂隙网络宏观非均质性；

裂隙内部开度微观尺度非均质性：基于逐次随机累加法，通过给定裂隙开度标准差σ和Hurst指数H，生成具有自仿射粗糙度特征的裂隙面，根据得到的每条裂隙的平均开度值b_correlated,i，使用位错法构建符合正态分布的随机裂隙开度场，并利用重心插值法将随机裂隙开度场投影到随机三维离散裂隙网络岩体模型每条裂隙的三角网格节点上，由此得到的裂隙内部开度差异即裂隙内部开度微观尺度非均质性；

步骤4.对裂隙面三角网格单元和基质四面体网格单元进行数据转换，得到具有多尺度非均质性且考虑岩石基质的三维裂隙网络岩体模型或网格模型。

进一步地，所述步骤1中蒙特卡洛模拟方法需要的参数为：

根据模拟对象实际情况确定的模拟区域的长、宽和高；

根据模拟对象实际情况确定的网络裂隙几何参数，包括数量、产状、密度和大小；

进一步地，步骤1中具体包括利用dfnWorks开源软件或其他建模软件或自主开发代码来使用蒙特卡洛模拟方法生成随机三维离散裂隙网络模型。

进一步地，所述步骤2需要将生成的裂隙几何参数进行处理转换得到Delaunay三角网格开源工具箱DFNsMeshGenerator3D的输入文件。

进一步地，裂隙开度与裂隙大小的幂律关系表达式为

b_correlated,i＝γr_i ^β

其中i为三维裂隙网络模型中第i条裂隙，i＝1，2，3，…，N，其中N为裂隙总数；b_correlated,i为第i条裂隙的平均裂隙开度；r_i为第i条裂隙的大小；γ和β为幂律关系的特征系数和特征指数。

进一步地，所述步骤3中构建裂隙内部开度微观尺度非均质性包括以下子步骤：

步骤3.3.1.基于逐次随机累加算法生成N个粗糙裂隙面，作为每条离散裂隙的下表面，其位置坐标用函数Z₁(i,X₀)表示，其中X₀＝(x,y)；

步骤3.3.2.复制下表面并根据需求剪切变形r_d进行位移生成上表面，其中r_d＝(x_d,y_d)，x_d和y_d分别表示x和y方向的剪切位移；生成的上表面位置坐标用函数Z₂(i,X₀)＝Z₁(i,X₀+r_d)+b_correlated,i表示

步骤3.3.3.得到每一条离散裂隙的二维笛卡尔坐标开度分布场，表示为

b_internal,i(X₀)＝Z₂(i,X₀)-Z₁(i,X₀)

利用MATLAB重心插值法在三维笛卡尔坐标中将生成的二维笛卡尔坐标裂隙开度场投影到三维裂隙网络模型每条裂隙的三角网格节点上；因此，三维裂隙网络模型中每条裂隙的裂隙开度分布b_internal,i(X)表示为

式中：X＝(x'，y'，z')，b_internal,i(X)＝0表示X坐标点的裂隙开度为0，即裂隙内部上下面接触区域。

进一步地，所述步骤3中构建裂隙内部开度微观尺度非均质性的H为Hurst指数，若裂隙为天然岩石，则取值范围为0.45～0.87，H值越小，粗糙度越大。

进一步地，所述步骤3中，在利用逐次随机累加算法算法生成二维笛卡尔坐标裂隙开度场时，需保证二维裂隙的中心与三维裂隙的中心对齐，且二维裂隙的尺寸大于三维裂隙网络模型中每一条裂隙的平面尺寸，以保证在使用MATLAB进行空间投影时开度场可完全覆盖裂隙网络模型中的裂隙；

进一步地，在利用MATLAB进行空间投影时，需根据每个网格结点坐标进行插值计算得到裂隙开度值，在此基础上对每个三角网格单元的3个结点取平均值，即可得到每个三角网格的平均裂隙开度，三角网格划分越精细，裂隙网络模型的裂隙开度场刻画越精确。

进一步地，所述步骤4中包括：

生成.vtk格式文件并使用Paraview对裂隙面三角网格单元和基质四面体网格单元进行数据转换，得到具有多尺度非均质性且考虑岩石基质的三维裂隙网络岩体模型或网格模型；根据需求将网格模型进行可视化分析，生成.msh格式文件并使用OpenGeosys开展渗流和溶质迁移模拟。

本发明的有益效果是：本发明不仅考虑了三维裂隙网络随机分布和复杂结构造成的宏观非均质性，而且还考虑了裂隙表面粗糙度引起的微观尺度非均质性，通过蒙特卡洛方法重构了与真实裂隙网络岩体相似的粗糙三维裂隙网络岩体模型，能够更加真实地模拟复杂裂隙网络结构及周围岩石基质，为岩体裂隙渗流和污染物迁移的理论研究工作提供了较为准确的数值模型。

附图说明

图1为本发明的流程图；

图2为本发明的裂隙网络结构示意图；

图3为本发明的裂隙面和边界面三角网格示意图；

图4为本发明的岩石基质四面体网格示意图；

图5为本发明的多尺度非均质性裂隙网络模型示意图；

图6为本发明的考虑多尺度非均质性和岩石基质的裂隙网络岩体模型示意图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

本发明提供一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，如图1所示，该方法包括以下步骤：

步骤1：首先，在Ubuntu系统下打开dfnWorks开源软件或其他建模软件并在输入文件中添加变量，设置模拟区域的长、宽和高均为15m，网络裂隙组数为1，裂隙数量为40，裂隙密度为0.5；裂隙产状呈Fisher分布，其中裂隙产状(倾向、倾角)和Fisher常数分别取值为90°、45°和10；裂隙大小呈截断幂函数分布，其中裂隙半径上下限取值分别为2m和8m，幂指数为1.8；裂隙平均开度与裂隙大小呈幂律关系，幂律关系的特征系数γ和特征指数β取值分别为5.0×10^-4和0.5。

其次，基于dfnWorks开源软件或其他建模软件的dfnGen模块，采用蒙特卡洛模拟方法生成如图2所示的随机三维离散裂隙网络模型。

步骤2：首先，利用Python脚本将dfnWorks生成的裂隙几何参数转换为一定格式的输入文件，并调用Delaunay三角网格开源工具箱DFNsMeshGenerator3D，同时对三维裂隙网络和基质边界面进行三角网格划分。为了保证在进行裂隙面三角网格划分时，同时获得6个边界面三角网格，在输入文件中需要增加6个裂隙面的顶点和中心坐标，将6个边界面作为裂隙面进行三角网格划分，裂隙面和边界面三角形网格划分结果如图3所示。图3中x、y、z为三维坐标系，黑色短实线为网格线段，三条短实线围合而成的三角形为上述三角网格。

然后，利用Python脚本提取裂隙面和基质边界面的三角网格单元构建PLC输入几何文件，包括三角网格单元的节点列表和面列表，并调用Delaunay四面体网格生成器TetGen对岩石基质进行四面体网格划分，如图4所示，图中为了展示四面体网格，将正方体裂隙岩体模型进行了切割处理，其中圆圈标示的区域展示了典型的四面体网格结构。

步骤3构建裂隙间宏观尺度非均质性、三维裂隙网络宏观尺度非均质性或裂隙内部开度微观尺度非均质性，并对随机三维离散裂隙网络岩体模型进行赋值或投影；构建非均质性包括以下子步骤：

构建裂隙间宏观尺度非均质性：首先对模拟区域中的所有裂隙进行编号，根据裂隙开度与裂隙大小之间的幂律关系，计算得到每条裂隙的平均开度值b_correlated,i；

裂隙开度与裂隙大小的幂律关系表达式为

b_correlated,i＝γr_i ^β

式中：i为随机三维离散裂隙网络模型中第i条裂隙，i＝1，2，3，…，N，其中N为裂隙总数，取值为40；b_correlated,i为第i条裂隙的裂隙开度；r_i为第i条裂隙的大小；γ和β为幂律关系的特征系数和特征指数。

对随机三维离散裂隙网络岩体模型每一条裂隙赋予每条裂隙的平均开度值b_correlated,i，由于每条裂隙大小不同，由此得到的裂隙间开度差异即裂隙间宏观尺度非均质性。

构建裂隙网络宏观尺度非均质性：根据随机三维离散裂隙网络模型中每一条裂隙的裂隙开度，可计算得到N条裂隙对应的平均开度值b_constant＝7.95e-5m，对随机三维离散裂隙网络模型每一条裂隙赋予N条裂隙对应的平均开度值b_constant，由此得到的随机三维离散裂隙网络模型可以用于模拟裂隙网络宏观非均质性。平均裂隙开度的计算表达式为

构建裂隙网络宏观尺度非均质性：首先，采用逐次随机累加法(SRA)生成具有自仿射粗糙度特征的下裂隙面，用函数Z₁(i,X₀)表示，X₀＝(x,y)。其中标准差(σ)取值为0.07mm，Hurst指数(H)取值为0.60。

其次，复制下裂隙面并通过剪切变形r_d＝(x_d,y_d)进行位移生成上表面，用函数Z₂(i,X₀)＝Z₁(i,X₀+r_d)+b_correlated,i表示，其中b_correlated,i为初始状态上下表面之间的平均裂隙开度。即根据随机三维离散裂隙网络模型中每条裂隙的平均裂隙开度先向上平移距离y_d＝b_correlated,i以模拟法向位移，再向右平移x_d＝1.5mm以模拟剪切作用下的水平位移，从而简单快速地构建符合正态分布的随机裂隙开度场。

然后，得到每一条离散裂隙的二维笛卡尔坐标开度分布场，可以表示为

b_internal,i(X₀)＝Z₂(i,X₀)-Z₁(i,X₀)。

利用MATLAB重心插值法在三维笛卡尔坐标中将生成的二维笛卡尔坐标裂隙开度场投影到随机三维离散裂隙网络模型每条裂隙的三角网格节点上，由此得到三维裂隙网络多尺度非均质性，如图5所示，图中展示了每条裂隙中的裂隙开度分布，Aperture含义为开度值，不同灰度代表不同开度值，且灰度越浅，开度值越小，灰度越深，开度值越大。其中随机三维离散裂隙网络模型中每条裂隙的裂隙开度分布b_internal,i(X)可以表示为

式中：X＝(x'，y'，z')，b_internal,i(X)＝0表示裂隙面内部的接触区域。

步骤4：利用Python脚本对裂隙面三角网格单元和基质四面体网格单元进行数据转换，可以得到具有多尺度非均质性且考虑岩石基质的三维裂隙网络岩体模型或网格模型，如图6所示，图中展示了具有不同裂隙开度的网络裂隙与周围四面体网格岩石基质组合而成的空间结构，Aperture含义为开度值。根据需求可以将网格模型转换为可视化文件或用于开展渗流和溶质迁移模拟的有限元网格，如生成.vtk格式文件并使用Paraview进行可视化分析，生成.msh格式文件并使用OpenGeosys(OGS)开展渗流和溶质迁移模拟。

本技术领域的人员根据本发明所提供的文字描述、附图以及权利要求书能够很容易在不脱离权利要求书所限定的本发明的思想和范围条件下，可以做出多种变化和改动。凡是依据本发明的技术思想和实质对上述实施例进行的任何修改、等同变化，均属于本发明的权利要求所限定的保护范围之内。

Claims (10)

1.一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，该方法包括以下步骤：

步骤1.利用蒙特卡洛模拟方法生成随机三维离散裂隙网络岩体模型；

步骤2.调用Delaunay三角网格开源工具箱DFNsMeshGenerator3D，对裂隙面和基质边界面进行三角网格划分，调用Delaunay四面体网格生成器TetGen对岩石基质进行四面体网格划分；

为了保证在进行裂隙面三角网格划分时，同时获得6个基质边界面的三角网格，在DFNsMeshGenerator3D输入文件中须增加6个裂隙面的顶点和中心坐标，将6个基质边界面进行三角网格划分；

为了保证在岩石基质四面体网格划分时，保持裂隙面及6个基质边界面的三角网格形态不变，TetGen的输入文件应包括裂隙面和基质边界面三角网格单元的节点列表和面列表；

步骤3.构建裂隙间宏观尺度非均质性、三维裂隙网络宏观尺度非均质性或裂隙内部开度微观尺度非均质性，并对随机三维离散裂隙网络岩体模型进行赋值或投影；构建非均质性包括以下子步骤：

构建裂隙间宏观尺度非均质性：对模拟区域中的所有裂隙进行编号，根据裂隙开度与裂隙大小之间的幂律关系，计算得到每条裂隙的平均开度值b_correlated,i，对随机三维离散裂隙网络岩体模型每一条裂隙赋予每条裂隙的平均开度值b_correlated,i，由于每条裂隙大小不同，由此得到的裂隙间开度差异即裂隙间宏观尺度非均质性；

构建裂隙网络宏观尺度非均质性：根据随机三维离散裂隙网络岩体模型中每条裂隙的平均裂隙开度b_correlated,i，计算得到所有裂隙的平均开度值b_constant，对随机三维离散裂隙网络岩体模型每一条裂隙的平均裂隙开度值赋予b_constant，由此得到的裂隙网络差异即裂隙网络宏观非均质性；

构建裂隙内部开度微观尺度非均质性：基于逐次随机累加法，通过给定裂隙开度标准差σ和Hurst指数H，生成具有自仿射粗糙度特征的裂隙面，根据得到的每条裂隙的平均开度值b_correlated,i，使用位错法构建符合正态分布的随机裂隙开度场，并利用重心插值法将随机裂隙开度场投影到随机三维离散裂隙网络岩体模型每条裂隙的三角网格节点上，由此得到的裂隙内部开度差异即裂隙内部开度微观尺度非均质性；

步骤4.对裂隙面三角网格单元和基质四面体网格单元进行数据转换，得到具有多尺度非均质性且考虑岩石基质的三维裂隙网络岩体模型或网格模型。

2.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤1中蒙特卡洛模拟方法需要的参数为：

根据模拟对象实际情况确定的模拟区域的长、宽和高；

根据模拟对象实际情况确定的网络裂隙几何参数，包括数量、产状、密度和大小。

3.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤1中具体包括利用dfnWorks开源软件或其他建模软件使用蒙特卡洛模拟方法生成随机三维离散裂隙网络模型。

4.根据权利要求2所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤2需要将生成的裂隙几何参数进行处理转换得到Delaunay三角网格开源工具箱DFNsMeshGenerator3D的输入文件。

5.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，裂隙开度与裂隙大小的幂律关系表达式为

b_correlated,i＝γr_i ^β

其中i为三维裂隙网络模型中第i条裂隙，i＝1，2，3，…，N，其中N为裂隙总数；b_correlated,i为第i条裂隙的平均裂隙开度；r_i为第i条裂隙的大小；γ和β为幂律关系的特征系数和特征指数。

6.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤3中构建裂隙内部开度微观尺度非均质性包括以下子步骤：

步骤3.3.1.基于逐次随机累加算法生成N个粗糙裂隙面，作为每条离散裂隙的下表面，其位置坐标用函数Z₁(i,X₀)表示，其中X₀＝(x,y)；

步骤3.3.2.复制下表面并根据需求剪切变形r_d进行位移生成上表面，其中r_d＝(x_d,y_d)，x_d和y_d分别表示x和y方向的剪切位移；生成的上表面位置坐标用函数Z₂(i,X₀)＝Z₁(i,X₀+r_d)+b_correlated,i表示；

步骤3.3.3.得到每一条离散裂隙的二维笛卡尔坐标开度分布场，表示为：

b_internal,i(X₀)＝Z₂(i,X₀)-Z₁(i,X₀)

利用MATLAB重心插值法在三维笛卡尔坐标中将生成的二维笛卡尔坐标裂隙开度场投影到三维裂隙网络模型每条裂隙的三角网格节点上；因此，三维裂隙网络模型中每条裂隙的裂隙开度分布b_internal,i(X)表示为：

式中：X＝(x'，y'，z')，b_internal,i(X)＝0表示X坐标点的裂隙开度为0，即裂隙内部上下面接触区域。

7.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤3中构建裂隙内部开度微观尺度非均质性的H为Hurst指数，若裂隙为天然岩石，则取值范围为0.45～0.87，H值越小，粗糙度越大。

8.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤3中，在利用逐次随机累加算法算法生成二维笛卡尔坐标裂隙开度场时，需保证二维裂隙的中心与三维裂隙的中心对齐，且二维裂隙的尺寸大于三维裂隙网络模型中每一条裂隙的平面尺寸，以保证在使用MATLAB进行空间投影时开度场可完全覆盖裂隙网络模型中的裂隙。

9.根据权利要求8所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，在利用MATLAB进行空间投影时，需根据每个网格结点坐标进行插值计算得到裂隙开度值，在此基础上对每个三角网格单元的3个结点取平均值，即可得到每个三角网格的平均裂隙开度，三角网格划分越精细，裂隙网络模型的裂隙开度场刻画越精确。

10.根据权利要求1所述的一种具有多尺度非均质性的三维裂隙网络岩体模型表征方法，其特征在于，所述步骤4中包括：

生成.vtk格式文件并使用Paraview对裂隙面三角网格单元和基质四面体网格单元进行数据转换，得到具有多尺度非均质性且考虑岩石基质的三维裂隙网络岩体模型或网格模型；根据需求将网格模型进行可视化分析，生成.msh格式文件并使用OpenGeosys开展渗流和溶质迁移模拟。

Images