CN117313589A - 三维粗糙离散裂隙网络构建、评价与渗流模拟的方法 - Google Patents
三维粗糙离散裂隙网络构建、评价与渗流模拟的方法 Download PDFInfo
- Publication number
- CN117313589A CN117313589A CN202311617961.2A CN202311617961A CN117313589A CN 117313589 A CN117313589 A CN 117313589A CN 202311617961 A CN202311617961 A CN 202311617961A CN 117313589 A CN117313589 A CN 117313589A
- Authority
- CN
- China
- Prior art keywords
- fracture
- rough
- fracture network
- model
- crack
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 70
- 238000009826 distribution Methods 0.000 claims abstract description 21
- 230000035699 permeability Effects 0.000 claims abstract description 11
- 238000004088 simulation Methods 0.000 claims abstract description 11
- 238000011835 investigation Methods 0.000 claims abstract description 8
- 238000005553 drilling Methods 0.000 claims abstract description 5
- 239000012530 fluid Substances 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 15
- 238000011161 development Methods 0.000 claims description 10
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 230000000694 effects Effects 0.000 claims description 4
- 238000004804 winding Methods 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 3
- 230000004069 differentiation Effects 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 238000005070 sampling Methods 0.000 claims description 2
- 239000011435 rock Substances 0.000 abstract description 13
- 238000004364 calculation method Methods 0.000 abstract description 8
- 238000011156 evaluation Methods 0.000 abstract description 7
- 238000011160 research Methods 0.000 description 11
- 229910000831 Steel Inorganic materials 0.000 description 4
- 239000010959 steel Substances 0.000 description 4
- 238000010276 construction Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000012800 visualization Methods 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
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/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- 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
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Abstract
本发明提供了一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，涉及裂隙岩体渗流计算技术领域，具体包括如下步骤：利用野外调查、钻孔数据、物探数据和其他实测数据获取裂隙网络几何参数概率分布模型。生成裂隙网络内各参数的数据集合。生成粗糙裂隙上、下壁面的各坐标下的高度值，并基于此获得所有裂隙的粗糙裂隙间开度集合。读取裂隙网络内各参数的数据集合，利用包围盒法进行三维随机裂隙网络连通性评价。建立三维离散、粗糙裂隙网络数值模型。计算三维裂隙网络数值模拟，求解模型渗透系数。本发明的技术方案克服现有技术中缺乏三维粗糙裂隙网络的建模方法、缺乏用于裂隙连通性评价与渗流模拟方法的问题。
Description
技术领域
本发明涉及裂隙岩体渗流计算技术领域，具体涉及一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法。
背景技术
天然岩体经历了漫长的成岩作用和复杂的地质构造演化，形成了大量的节理裂隙，显著改变了岩体的力学性质和渗流特性，影响着地下工程的建设安全，因此，基于野外裂隙产状建立真实的三维离散-粗糙裂隙网络是开展地下工程灾害防控研究的重点。
近年来，随着计算机技术和数值计算方法的发展，数值计算越来越多地应用于岩石力学岩石工程问题中，从岩体结构与形态出发，基于大量的野外岩体裂隙实测资料统计，通过概率论和数理统计理论进行离散裂隙网络建模(DFN)的方法在众多实际工程中得到了应用与验证，二维离散裂隙网络模型生成、可视化以及渗透系数求解等方法与流程被广泛研究。然而，真实世界中，裂隙网络以三维裂隙存在，二维离散裂隙网络模型存在误差，为此，Priest、王恩志、张国强、谢静等开展了三维裂隙网络的建模，通过生成裂隙参数（密度、长度、倾向、倾角、裂隙开度）借助圆盘模型、四边形模型生成三维裂隙网络。同时，关于三维裂隙网络渗流特性的数值计算研究得以发展。一方面，学者们借助基于CFD方法的主流数值模拟软件，开展渗流特性研究，例如Ansys Fluent，OpenFOAM，Comsol Multiphysics等，该类商业软件开发程度高，理论完备，基本可以满足大多数的模拟需求；另一方面，借助C++，Java，Matlab等编程平台，自行开发可用的三维渗流模拟软件开展研究，当前比较流行的方法有伽辽金法、3D Unified pipe-network Method(UPM)以及VOF法。
然而当前研究存在三个方面问题:(1) 三维裂隙多假设被光滑平板裂隙，然而真实裂隙是粗糙的，裂隙的粗糙程度将显著改变单一、网络裂隙的渗透性，现有研究无法考虑裂隙面的粗糙性；（2）三维裂隙岩体中裂隙随机分布，往往具有分布变异性大、数量多的特点，相较两种裂隙模拟方法，当裂隙数量过多时，商业软件往往无法对其直接建模，在裂隙属性赋值、边界条件设置上存在众多困难，需进行必要的二次开发；自行开发的仿真程序计算速度快，收敛性强，但往往存在开发程度低，应用范围小等难题。（3）裂隙网络渗流特性过程的研究包括，裂隙网络的生成、连通性的计算以及渗透系数的求解，面对复杂的真实裂隙网络世界，一次求解单个裂隙网络模型效率低，缺少可进行大规模自动求解的方法。
因此，现需要一种基于岩层产状的三维离散粗糙裂隙岩体建模、裂隙连通性评价、渗流模拟的较完整的技术方法，以便于开展三维粗糙裂隙网络的建模、评价与模拟，服务实际水利工程。
发明内容
本发明的主要目的在于提供一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，以解决现有技术中缺乏三维粗糙裂隙网络的建模方法、缺乏用于裂隙连通性评价与渗流模拟方法的问题。
为实现上述目的，本发明提供了一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，具体包括如下步骤：
S1，利用野外调查、钻孔数据、物探数据和其他实测数据获取裂隙网络几何参数概率分布模型。
S2，利用matlab编程软件，基于Monte-Carlo法对裂隙圆盘圆心、裂隙倾
角、裂隙倾向、裂隙圆盘半径和裂隙开度参数进行
随机抽样，生成裂隙网络内各参数的数据集合。
S3，使用随机形式的维尔斯特拉斯函数来生成粗糙裂隙上、下壁面的各坐标下的高度值，并基于此获得所有裂隙的粗糙裂隙间开度集合。
S4，读取由步骤S2获取的裂隙网络内各参数的数据集合，利用包围盒法进行三维随机裂隙网络连通性评价。
S5，基于Matlab with Comsol开发平台，利用欧拉角与旋转矩阵方法进行裂隙圆盘的坐标系转换与粗糙裂隙开度赋值，建立三维离散、粗糙裂隙网络数值模型。
S6，对三维裂隙网络渗流进行模拟，并计算三维裂隙网络模型渗透系数。
进一步地，步骤S1中裂隙网络几何参数分布模型包括：裂隙圆盘圆心模型、裂隙倾角模型、裂隙倾向模型、裂隙圆盘半径模型和裂隙开度模型。
进一步地，步骤S2中裂隙网络内各个数据集合，具体为：
；
；
；
；
；
。
式中：为裂隙模型内的裂隙数量，、和分别代表裂隙模型的长、宽、高；代表第个裂隙中心点在笛卡尔坐标系中的坐标，为第个裂隙的倾角，
为第个裂隙的倾向，为第个裂隙圆盘的半径，为第个裂隙的开度，为裂隙模型
裂隙密度。
进一步地，步骤S3具体包括如下步骤：
S3.1，使用随机形式的维尔斯特拉斯函数来生成粗糙裂隙上、下壁面的各坐标下的高度值，其中上裂隙壁面高度为:
（1）。
式中，为上粗糙裂隙壁面内中的平面坐标高度，是
具有标准正态分布的独立随机数，是随机数的个数，和是分形变量，和是具有
[0，2π]均匀分布的独立随机数，, 为生成裂隙网络的空间频率分辨率；改变的随机数
种子，以生成具有差异的粗糙裂隙上、下壁面数据集与 。
S3.2，读取步骤S2中裂隙网络内的裂隙开度集合，
令, 即将生成的裂隙数据集
赋予初始开度为的裂隙，为未赋予初始开度的裂隙数据集。
S3.3，考虑地应力作用，通过输入应力、法向刚度获取裂隙闭合变形量，获取实际粗糙裂隙开度变形量与上、下壁面接触情况：
（2）；
式中，为裂隙闭合变形量，为初始裂隙闭合变形，是正应力，是初始
法相刚度。
S3.4，上、下裂隙位移部分闭合后的粗糙裂隙平均裂隙开度利用公式（3）计算：
（3）；
其中，和分别代表x方向和y方向的节点的个数。
S3.5，循环步骤S3.1~ S3.4，以获得所有裂隙的粗糙裂隙间开度集合。
进一步地，步骤S4具体包括如下步骤：
S4.1，利用包围盒法对第一裂隙和第二裂隙的空间关系进行初步判断，对比两包
围球体球心间的直线距离与球体半径分析两裂隙包围球体重叠情况，几何特征满足公
式（4）：
（4）。
式中：为第一裂隙和第二裂隙包围球之间的距离，（，，）、，
，分别为第一裂隙圆盘、第二裂隙圆盘的圆心全局坐标，和分别表示第一裂隙圆
盘、第二裂隙圆盘的半径。
S4.2，当时，裂隙圆盘法向量表示为：
（5）。
为方向向量；为倾向，为倾角。
S4.3，第一个裂隙圆盘所在的平面表示为：
（6）；
其中，、、为第一个裂隙圆盘的方向向量。
第一个裂隙圆盘的边界轮廓线表示为：
（7）。
同理，第二个裂隙圆盘所在平面及边界轮廓线可分别表示为：
（8）。
其中，、、为第二个裂隙圆盘的方向向量。
（9）。
S4.4，当两个包围球体部分重叠后，利用公式（10）计算两裂隙圆盘所在平面的夹
角：
（10）。
S4.5，当夹角°时，表明两裂隙圆盘所在平面平行，裂隙圆盘交线不存在；相
反地，则分别计算两个裂隙圆盘是否与交线相交。
进一步地，步骤S5具体包括如下步骤：
S5.1，在Matlab中读取加载步骤S2中裂隙网络几何参数集合，以平面全局坐标系
原点为裂隙圆盘圆心，建立半径为的裂隙圆盘。
S5.2，读取步骤S3中新生成的上、下粗糙裂隙面间的裂隙开度值，在
Matlab with Comsol开发平台中输入命令，将其赋值于裂隙圆盘中，使得裂隙圆盘具有粗
糙裂隙开度。
S5.3，将全局坐标系的原点平移到圆盘中心点，即将全局坐标系原点沿三个方向的平移方向分别为。
S5.4，引入欧拉角定理，将坐标系绕轴旋转，使得全局坐标系的轴
旋转至裂隙的倾向方向；将坐标系绕轴旋转，使得全局坐标系的轴与裂隙的倾向线重
合，在笛卡尔坐标系中建立裂隙圆盘模型。
S5.5，重复步骤S5.1~ S5.4，逐次迭代，完成三维空间内的离散、粗糙裂隙网络建模。
进一步地，步骤S6具体包括如下步骤：
S6.1，基于matlab with comsol软件平台，调用软件平台的控制方程与求解器，其中在三维裂隙网络渗流模拟过程中的稳态流动由Darcy方程控制，具体为：
（11）；
（12）；
（13）；
式中：为裂隙开度，/>为裂隙密度，/>为动力黏度，/>为裂隙流量，/>为速度，/>为压力，/>为渗透率；/>为微积分的梯度算子，代表在不同方向上的微分。
S6.2，设定边界条件，以方向作为流体流动方向，在模型入口与出口边界分别设
置为压力边界条件。
S6.3，利用comsol模拟软件进行模型网格划分，并对流场和压力场进行求解。
S6.4，计算裂隙网络模型入口渗流速度，通过在入口处对流速积分，获取不同裂隙
开度分布情况下的稳态体积流速，具体为：
（14）。
式中：，为垂直于流体流动方向的横截面积，为流体通过横截面积的
速度。
S6.5，根据达西方程，计算模型的裂隙网络模型渗透系数，具体为：
（15）。
式中：为垂直于流体流动方向的横截面积，为裂隙网络在流动方向上的长度，为模型进出口边界压差。
本发明具有如下有益效果：
基于岩石裂隙产状生成裂隙网络参数，借助欧拉角与旋转矩阵解决了comsol商业软件平台内的裂隙圆盘的坐标系转换问题。
实现了三维“离散-粗糙”裂隙网络的建模与可视化。
一次性自动实现从裂隙网络建模到裂隙网络连通性评价，再到渗流模拟的过程，降低人工操作成本，提高计算效率，使得大规模数值模拟研究成为可能。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1示出了本发明的三维粗糙离散裂隙网络构建、评价与渗流模拟的方法的流程图。
图2示两个包围球为相离状态的示意图。
图3示出了两个包围球为相交状态示意图。
图4示出了裂隙相互包含的状态示意图。
图5示出了裂隙相交的状态示意图。
图6示出了裂隙相离的状态示意图。
图7示出了根据步骤S5.1建立的三维离散裂隙网络模型图。
图8示出了根据步骤S5.2得到的粗糙裂隙的裂隙开度赋值示意图。
具体实施方式
下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，具体包括如下步骤：
S1，利用野外调查、钻孔数据、物探数据和其他实测数据获取裂隙网络几何参数概率分布模型，如表1所示。
表1：裂隙网络几何参数概率分布模型
。
S2，利用matlab编程软件，基于Monte-Carlo法对裂隙圆盘圆心、裂隙倾
角、裂隙倾向、裂隙圆盘半径和裂隙开度参数进行
随机抽样，生成裂隙网络内各参数的数据集合，如表2所示。
表2：裂隙网络内各参数的数据集合
。
S3，使用随机形式的维尔斯特拉斯函数来生成粗糙裂隙上、下壁面的各坐标下的高度值，并基于此获得所有裂隙的粗糙裂隙间开度集合。
S4，读取由步骤S2获取的裂隙网络内各参数的数据集合，利用包围盒法进行三维随机裂隙网络连通性评价。
S5，基于Matlab with Comsol开发平台，利用欧拉角与旋转矩阵方法进行裂隙圆盘的坐标系转换与粗糙裂隙开度赋值，建立三维离散、粗糙裂隙网络数值模型。
S6，对三维裂隙网络渗流进行模拟，并计算三维裂隙网络模型渗透系数。
具体地，开展现场勘察调研，利用野外调查数据、钻孔数据、物探数据等实测数据获取裂隙网络参数概率分布模型，确定各概率分布模型的参数值。本发明采用Baecher裂隙圆盘模型，步骤S1中裂隙网络几何参数分布模型包括：裂隙圆盘圆心模型、裂隙倾角模型、裂隙倾向模型、裂隙圆盘半径模型和裂隙开度模型。
裂隙圆盘圆心模型符合泊松分布，公式为：
。
裂隙倾角模型、裂隙倾向模型符合Fisher分布，公式为：
。
裂隙圆盘半径模型符合幂律分布，公式为：
。
裂隙开度模型符合均匀分布函数，公式为：
。
具体地，步骤S2中裂隙网络内各个数据集合，具体为：
；
；
；
；
；
。
式中：为裂隙模型内的裂隙数量，、和分别代表裂隙模型的长、宽、高，代表第个裂隙中心点在笛卡尔坐标系中的坐标，为第个裂隙的倾角，
为第个裂隙的倾向，为第个裂隙圆盘的半径，为第个裂隙的开度，为裂隙模型
裂隙密度。
具体地，现实世界中裂隙是粗糙的，裂隙上、下粗糙壁面间的裂隙空间（即裂隙开度）影响着裂隙内的渗流特性。目前大多数的粗糙裂隙研究多是针对单一裂隙研开展的，三维离散裂隙网络内粗糙裂隙建模与渗流特征研究相对较少。生成三维离散粗糙裂隙网络，首先要为每个裂隙生成各自的粗糙裂隙间开度数据集合，以单一粗糙裂隙生成过程为例，步骤S3具体包括如下步骤：
S3.1，为了使得生成的裂隙开度更符合实际，需先生成单一裂隙的上、下粗糙裂隙壁面。裂隙面的粗糙程度可由分形维数表征，使用随机形式的维尔斯特拉斯函数来生成粗糙裂隙上、下壁面的各坐标下的高度值，其中上裂隙壁面高度为:
（1）；
式中，为上粗糙裂隙壁面内中的平面坐标高度，是
具有标准正态分布的独立随机数，是随机数的个数，和是分形变量，和是具有
[0，2π]均匀分布的独立随机数，, 为生成裂隙网络的空间频率分辨率；通常取64、128、
256、512、1024等。
在粗糙裂隙上、下壁面生成过程中，可在其他参数保持不变的情况下，改变的
随机数种子，以生成具有差异的粗糙裂隙上、下壁面数据集与。
S3.2，读取步骤S2中裂隙网络内的裂隙开度集合，
令, 即将生成的裂隙数据集
赋予初始开度为的裂隙，为未赋予初始开度的裂隙数据集。
S3.3，考虑地应力作用，通过输入应力、法向刚度获取裂隙闭合变形量，获取实际粗糙裂隙开度变形量与上、下壁面接触情况：
（2）；
式中，为裂隙闭合变形量，为初始裂隙闭合变形，是正应力，是初始
法相刚度。
S3.4，令降低（下移）。即=
-；重新计算,即-，当存在时，说明此时上、下裂隙壁面接触，
该位置裂隙闭合，流体无法通过，则令=0。更新数据，即可获得单一粗糙裂隙间开
度数据集合。上、下裂隙位移，部分闭合后的粗糙裂隙平均裂隙开度利用公式（3）计算：
（3）；
其中，和分别代表x方向和y方向的节点的个数。
S3.5，循环步骤S3.1~ S3.4，以获得所有裂隙的粗糙裂隙间开度集合。
具体地，读取各裂隙网络参数集合，进行三维随机裂隙网络连通性评价。由于在裂
隙岩体渗流过程中，孤立裂隙无法对流体的运移过程起到实质作用，因此流体在裂隙网络
中的流动主要取决于岩体裂隙网络的连通性。裂隙连通性表征为裂隙交叉点、线的个数，连
通性指数可定义为：
式中：研究区为三维对象时，为区域内裂隙圆盘间相交线段的个数，为区域内
裂隙圆盘的个数。
步骤S4具体包括如下步骤：
S4.1，为计算裂隙圆盘相交频次，利用包围盒法对第一裂隙和第二裂隙的空间关
系进行初步判断，对比两包围球体球心间的直线距离与球体半径分析两裂隙包围球体
重叠情况，如图2和图3所示，几何特征满足公式（4）：
（4）。
式中：为第一裂隙和第二裂隙包围球之间的距离，（，，）、，
，分别为第一裂隙圆盘、第二裂隙圆盘的圆心全局坐标，和分别表示第一裂隙圆
盘、第二裂隙圆盘的半径。
S4.2，当时，应进一步对两裂隙位置关系进行计算判断。根据裂
隙圆盘产状（倾角和倾向，可根据倾角弧度值、倾向弧值值转换）裂隙圆盘法向量表示为：
（5）；
为方向向量；为倾向，为倾角。
S4.3，第一个裂隙圆盘所在的平面表示为：
（6）；
其中，、、为第一个裂隙圆盘的方向向量。
第一个裂隙圆盘的边界轮廓线表示为：
（7）。
同理，第二个裂隙圆盘所在平面及边界轮廓线可分别表示为：
（8）。
其中，、、为第二个裂隙圆盘的方向向量；
（9）。
S4.4，当两个包围球体部分重叠后，利用公式（10）计算两裂隙圆盘所在平面的夹
角：
（10）。
S4.5，当夹角时，表明两裂隙圆盘所在平面平行，裂隙圆盘交线不存在；相
反地，则利用公式（5）-公式（10）分别计算两个裂隙圆盘是否与交线相交。当两个裂隙圆盘
均与交线相交时，应利用两裂隙圆盘与交线的四个空间交点，分析判断裂隙的位置关系。如
图4和图5所示，第二个裂隙被包含于或相交于第一个裂隙，说明两个裂隙圆盘存在相交关
系，与此同时，可利用四个坐标求得交线空间方程及相交长度；图6中两裂隙相离，表明两个
裂隙圆盘虽有统一的交线，但两裂隙圆盘不相交。其中A11 、A12 分别为裂隙A与交线的第一
个、第二个交点；B11 和B12分别为裂隙B与交线的第一个、第二个交点。
具体地，以第一个裂隙圆盘生成过程为例，假设第一个裂隙圆盘参数分别为，倾角与倾向，圆盘半径。
步骤S5具体包括如下步骤：
S5.1，在Matlab中读取加载步骤S2中裂隙网络几何参数集合，以平面全局坐标系
原点为裂隙圆盘圆心，建立半径为的裂隙圆盘，如图7所示实现三维离散裂隙网络模型的
可视化。
S5.2，读取步骤S3中新生成的上、下粗糙裂隙面间的裂隙开度值，在
Matlab with Comsol开发平台中输入命令，将其赋值于裂隙圆盘中，使得裂隙圆盘具有粗
糙裂隙开度，如图8所示，即可实现粗糙裂隙的裂隙开度赋值与可视化。
S5.3，将全局坐标系的原点平移到圆盘中心点，即将全局坐标系原点沿三个方向的平移方向分别为。
S5.4，引入欧拉角定理，将坐标系绕轴旋转，使得全局坐标系的轴
旋转至裂隙的倾向方向；将坐标系绕轴旋转，使得全局坐标系的轴与裂隙的倾向线重
合，在笛卡尔坐标系中建立裂隙圆盘模型。
S5.5，重复步骤S5.1~ S5.4，逐次迭代，完成三维空间内的离散、粗糙裂隙网络建模。
具体地，步骤S6具体包括如下步骤：
S6.1，基于matlab with comsol软件平台，调用软件平台的控制方程与求解器，其中在三维裂隙网络渗流模拟过程中的稳态流动由Darcy方程控制，具体为：
（11）；
（12）；
（13）；
式中：为裂隙开度，/>为裂隙密度，/>为动力黏度，/>为裂隙流量，/>为速度，/>为压力，/>为渗透率；/>为微积分的梯度算子，代表在不同方向上的微分。
S6.2，设定边界条件，以方向作为流体流动方向，在模型入口与出口边界分别设
置为压力边界条件。
S6.3，利用comsol模拟软件进行模型网格划分，并对流场和压力场进行求解。
S6.4，计算裂隙网络模型入口渗流速度，通过在入口处对流速积分，获取不同裂隙
开度分布情况下的稳态体积流速，具体为： （14）；
式中：，为垂直于流体流动方向的横截面积，为流体通过横截面积的
速度。
S6.5，根据达西方程，计算模型的裂隙网络模型渗透系数，具体为：
（15）；
式中：为垂直于流体流动方向的横截面积，为裂隙网络在流动方向上的长度，为模型进出口边界压差。
以上为单个三维粗糙裂隙网络建模、连通性评价以及渗流求解过程。通过编写脚本程序，不断读取试验方案，生成裂隙网络参数集合，进行程序循环。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
Claims (7)
1.一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，具体包括如下步骤：
S1，利用野外调查、钻孔数据、物探数据和其他实测数据获取裂隙网络几何参数概率分布模型；
S2，利用matlab编程软件，基于Monte-Carlo法对裂隙圆盘圆心、裂隙倾角、裂隙倾向/>、裂隙圆盘半径/>和裂隙开度/>参数进行随机抽样，生成裂隙网络内各参数的数据集合；
S3，使用随机形式的维尔斯特拉斯函数来生成粗糙裂隙上、下壁面的各坐标下的高度值，并基于此获得所有裂隙的粗糙裂隙间开度集合；
S4，读取由步骤S2获取的裂隙网络内各参数的数据集合，利用包围盒法进行三维随机裂隙网络连通性评价；
S5，基于Matlab with Comsol开发平台，利用欧拉角与旋转矩阵方法进行裂隙圆盘的坐标系转换与粗糙裂隙开度赋值，建立三维离散、粗糙裂隙网络数值模型；
S6，对三维裂隙网络渗流进行模拟，并计算三维裂隙网络模型渗透系数。
2.根据权利要求1所述的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，步骤S1中裂隙网络几何参数分布模型包括：裂隙圆盘圆心模型、裂隙倾角模型、裂隙倾向模型、裂隙圆盘半径模型和裂隙开度模型。
3.根据权利要求1所述的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，步骤S2中裂隙网络内各个数据集合，具体为：
；
；
；
；
；
；
式中：为裂隙模型内的裂隙数量，/>、/>和/>分别代表裂隙模型的长、宽、高；代表第/>个裂隙中心点在笛卡尔坐标系中的坐标，/>为第/>个裂隙的倾角，/>为第/>个裂隙的倾向，/>为第/>个裂隙圆盘的半径，/>为第/>个裂隙的开度，/>为裂隙模型裂隙密度。
4.根据权利要求1所述的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，步骤S3具体包括如下步骤：
S3.1，使用随机形式的维尔斯特拉斯函数来生成粗糙裂隙上、下壁面的各坐标下的高度值，其中上裂隙壁面高度为:
（1）；
式中，为上粗糙裂隙壁面内/>中的平面坐标高度，/>是具有标准正态分布的独立随机数，/>是随机数的个数，/>和/>是分形变量，/>和/>是具有[0，2π]均匀分布的独立随机数，/>, />为生成裂隙网络的空间频率分辨率；改变/>的随机数种子，以生成具有差异的粗糙裂隙上、下壁面数据集/>与/> ；
S3.2，读取步骤S2中裂隙网络内的裂隙开度集合，令/>, 即将生成的裂隙数据集/>赋予初始开度为/>的裂隙，/>为未赋予初始开度的裂隙数据集；
S3.3，考虑地应力作用，通过输入应力、法向刚度获取裂隙闭合变形量，获取实际粗糙裂隙开度变形量与上、下壁面接触情况：
（2）；
式中，为裂隙闭合变形量，/>为初始裂隙闭合变形，/>是正应力，/>是初始法相刚度；
S3.4，上、下裂隙位移部分闭合后的粗糙裂隙平均裂隙开度利用公式（3）计算：
（3）；
其中，和/>分别代表x方向和y方向的节点的个数；
S3.5，循环步骤S3.1~ S3.4，以获得所有裂隙的粗糙裂隙间开度集合。
5.根据权利要求1所述的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，步骤S4具体包括如下步骤：
S4.1，利用包围盒法对第一裂隙和第二裂隙的空间关系进行初步判断，对比两包围球体球心间的直线距离与球体半径分析两裂隙包围球体重叠情况，几何特征满足公式（4）：
（4）；
式中：为第一裂隙和第二裂隙包围球之间的距离，/>（/>，/>，/>）、/>，/>，/>分别为第一裂隙圆盘、第二裂隙圆盘的圆心全局坐标，/>和/>分别表示第一裂隙圆盘、第二裂隙圆盘的半径；
S4.2，当时，裂隙圆盘法向量/>表示为：
（5）；
为方向向量；/>为倾向，/>为倾角；
S4.3，第一个裂隙圆盘所在的平面表示为：
（6）；
其中，、/>、/>为第一个裂隙圆盘的方向向量；
第一个裂隙圆盘的边界轮廓线表示为：
（7）；
同理，第二个裂隙圆盘所在平面及边界轮廓线可分别表示为：
（8）；
其中，、/>、/>为第二个裂隙圆盘的方向向量；
（9）；
S4.4，当两个包围球体部分重叠后，利用公式（10）计算两裂隙圆盘所在平面的夹角：
（10）；
S4.5，当夹角时，表明两裂隙圆盘所在平面平行，裂隙圆盘交线不存在；相反地，则分别计算两个裂隙圆盘是否与交线相交。
6.根据权利要求1所述的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，步骤S5具体包括如下步骤：
S5.1，在Matlab中读取加载步骤S2中裂隙网络几何参数集合，以平面全局坐标系原点为裂隙圆盘圆心，建立半径为的裂隙圆盘；
S5.2，读取步骤S3中新生成的上、下粗糙裂隙面间的裂隙开度值，在Matlabwith Comsol开发平台中输入命令，将其赋值于裂隙圆盘中，使得裂隙圆盘具有粗糙裂隙开度；
S5.3，将全局坐标系的原点平移到圆盘中心点，即将全局坐标系原点沿/>三个方向的平移方向分别为/>；
S5.4，引入欧拉角定理，将坐标系绕轴旋转/>，使得全局坐标系的/>轴旋转至裂隙的倾向方向；将坐标系绕/>轴旋转/>，使得全局坐标系的/>轴与裂隙的倾向线重合，在笛卡尔坐标系中建立裂隙圆盘模型；
S5.5，重复步骤S5.1~ S5.4，逐次迭代，完成三维空间内的离散、粗糙裂隙网络建模。
7.根据权利要求1所述的一种三维粗糙离散裂隙网络构建、评价与渗流模拟的方法，其特征在于，步骤S6具体包括如下步骤：
S6.1，基于matlab with comsol软件平台，调用软件平台的控制方程与求解器，其中在三维裂隙网络渗流模拟过程中的稳态流动由Darcy方程控制，具体为：
（11）；
（12）；
（13）；
式中：为裂隙开度，/>为裂隙密度，/>为动力黏度，/>为裂隙流量，/>为速度，/>为压力，/>为渗透率；/>为微积分的梯度算子，代表在不同方向上的微分；
S6.2，设定边界条件，以方向作为流体流动方向，在模型入口与出口边界分别设置为压力边界条件；
S6.3，利用comsol模拟软件进行模型网格划分，并对流场和压力场进行求解；
S6.4，计算裂隙网络模型入口渗流速度，通过在入口处对流速积分，获取不同裂隙开度分布情况下的稳态体积流速，具体为：
（14）；
式中：，为垂直于流体流动方向的横截面积，/>为流体通过横截面积的速度；
S6.5，根据达西方程，计算模型的裂隙网络模型渗透系数，具体为：
（15）；
式中：为垂直于流体流动方向的横截面积，/>为裂隙网络在流动方向上的长度，/>为模型进出口边界压差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311617961.2A CN117313589B (zh) | 2023-11-30 | 2023-11-30 | 三维粗糙离散裂隙网络构建、评价与渗流模拟的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311617961.2A CN117313589B (zh) | 2023-11-30 | 2023-11-30 | 三维粗糙离散裂隙网络构建、评价与渗流模拟的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117313589A true CN117313589A (zh) | 2023-12-29 |
| CN117313589B CN117313589B (zh) | 2024-02-06 |
Family
ID=89281603
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311617961.2A Active CN117313589B (zh) | 2023-11-30 | 2023-11-30 | 三维粗糙离散裂隙网络构建、评价与渗流模拟的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117313589B (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106991244A (zh) * | 2017-04-13 | 2017-07-28 | 河海大学 | 一种基于图论的裂隙网络连通性及渗流计算的方法 |
| CN112131642A (zh) * | 2020-09-17 | 2020-12-25 | 四川大学 | 一种基于波数法的岩体粗糙离散裂隙网络生成方法 |
| AU2020103292A4 (en) * | 2020-11-06 | 2021-01-14 | Karimzade, Ebrahim DR | A novel approach for evaluation of the flow paths into underground excavations using three-dimensional discrete fracture network (3D DFN) model |
| WO2021022465A1 (zh) * | 2019-08-02 | 2021-02-11 | 山东科技大学 | 一种基于数字图像技术的粗糙面裂隙生成方法及试验系统 |
| WO2022011894A1 (zh) * | 2020-07-15 | 2022-01-20 | 中海油田服务股份有限公司 | 基于卷积神经网络的孔隙网络模型的建模方法及装置 |
| US20220136387A1 (en) * | 2020-11-02 | 2022-05-05 | Southwest Petroleum University | Dynamic fracture width calculation method for drilling fluid loss in fractured formation |
| CN115310381A (zh) * | 2022-08-08 | 2022-11-08 | 中国电建集团成都勘测设计研究院有限公司 | 基于典型单元体计算的裂隙岩体渗流场模拟方法及终端 |
| CN115661388A (zh) * | 2022-11-02 | 2023-01-31 | 浙大城市学院 | 具有多尺度非均质性的三维裂隙网络岩体模型表征方法 |
| CN116630243A (zh) * | 2023-04-24 | 2023-08-22 | 中建三局第一建设工程有限责任公司 | 裂隙岩体渗流过程模拟方法、装置、计算机设备及介质 |
-
2023
- 2023-11-30 CN CN202311617961.2A patent/CN117313589B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106991244A (zh) * | 2017-04-13 | 2017-07-28 | 河海大学 | 一种基于图论的裂隙网络连通性及渗流计算的方法 |
| WO2021022465A1 (zh) * | 2019-08-02 | 2021-02-11 | 山东科技大学 | 一种基于数字图像技术的粗糙面裂隙生成方法及试验系统 |
| WO2022011894A1 (zh) * | 2020-07-15 | 2022-01-20 | 中海油田服务股份有限公司 | 基于卷积神经网络的孔隙网络模型的建模方法及装置 |
| CN112131642A (zh) * | 2020-09-17 | 2020-12-25 | 四川大学 | 一种基于波数法的岩体粗糙离散裂隙网络生成方法 |
| US20220136387A1 (en) * | 2020-11-02 | 2022-05-05 | Southwest Petroleum University | Dynamic fracture width calculation method for drilling fluid loss in fractured formation |
| AU2020103292A4 (en) * | 2020-11-06 | 2021-01-14 | Karimzade, Ebrahim DR | A novel approach for evaluation of the flow paths into underground excavations using three-dimensional discrete fracture network (3D DFN) model |
| CN115310381A (zh) * | 2022-08-08 | 2022-11-08 | 中国电建集团成都勘测设计研究院有限公司 | 基于典型单元体计算的裂隙岩体渗流场模拟方法及终端 |
| CN115661388A (zh) * | 2022-11-02 | 2023-01-31 | 浙大城市学院 | 具有多尺度非均质性的三维裂隙网络岩体模型表征方法 |
| CN116630243A (zh) * | 2023-04-24 | 2023-08-22 | 中建三局第一建设工程有限责任公司 | 裂隙岩体渗流过程模拟方法、装置、计算机设备及介质 |
Non-Patent Citations (3)
| Title |
|---|
| 何忱;姚池;杨建华;蒋水华;姜清辉;周创兵;: "基于等效离散裂隙网络的三维裂隙岩体渗流模型", 岩石力学与工程学报, no. 1 * |
| 宋晓晨, 徐卫亚: "裂隙岩体渗流模拟的三维离散裂隙网络数值模型(Ⅰ):裂隙网络的随机生成", 岩石力学与工程学报, no. 12 * |
| 张亚;彭乐文;: "基于Monte-Carlo的岩体三维裂隙网络模型研究", 甘肃科学学报, no. 02 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117313589B (zh) | 2024-02-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107060746B (zh) | 一种复杂裂缝性油藏流动模拟的方法 | |
| CN102254349B (zh) | 一种使用钻孔数据构建沉积地层系统三维实体模型的方法 | |
| CN109241588B (zh) | 一种基于拟连续地质力学模型的单裂缝扩展的模拟方法 | |
| CN106960070A (zh) | 一种基于有限元‑离散元ct重构煤体的渗流模拟方法 | |
| CN105225273B (zh) | 基于空间断层多边形建立断层模型的方法及装置 | |
| CN107203670B (zh) | 一种粗糙随机型离散节理网络模型构建方法 | |
| CN106934826B (zh) | 一种岩质边坡结构精细化建模及块体识别方法 | |
| CN102640155A (zh) | 计算储层模拟器的近似油井泄油压力的系统、计算机实现方法和计算机可读程序产品 | |
| Xue et al. | A fast numerical method and optimization of 3D discrete fracture network considering fracture aperture heterogeneity | |
| CN113189305B (zh) | 一种孔隙型可渗透岩石的仿真模拟方法及系统 | |
| CN112326524A (zh) | 一种基于ct扫描图像的岩石孔渗测量方法 | |
| CN112069654B (zh) | 一种碳酸盐岩酸化数值模拟方法 | |
| CN114549774A (zh) | 一种基于钻孔数据的三维地层建模方法 | |
| Zhao et al. | Three-dimensional representation of discrete fracture matrix model for fractured reservoirs | |
| Zeng et al. | Refined modeling and movement characteristics analyses of irregularly shaped particles | |
| CN115238553A (zh) | 一种地埋管线渗漏浸蚀的危险区域划分方法和系统 | |
| Chen et al. | Spatial estimation of material parameters and refined finite-element analysis of rockfill dam based on construction digitization | |
| CN117313589B (zh) | 三维粗糙离散裂隙网络构建、评价与渗流模拟的方法 | |
| Ma et al. | Study on representative volume elements considering inhomogeneity and anisotropy of rock masses characterised by non-persistent fractures | |
| EP2678718A2 (en) | Computerized method for the estimation of a value for at least a parameter of a hydrocarbon-producing region, for planning the operation and operating the region | |
| CN110990926B (zh) | 一种基于面积修正率的城市地表建筑水动力学仿真方法 | |
| CN107507179A (zh) | 基于gocad的岩土体量化分析方法 | |
| CN109472046A (zh) | 复杂坝基拱坝三维有限元四面体网格自动剖分方法 | |
| CN108710769A (zh) | 一种离散元仿真中不规则砂土颗粒模型的构建方法 | |
| CN113530536A (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 |