一种基于图论的裂隙网络连通性及渗流计算的方法与流程

本发明涉及一种模拟裂隙连通性及渗流的方法，具体是一种基于图论的裂隙网络连通性及渗流计算的方法。

背景技术：

断裂的地质体在世界范围内普遍存在，人们很早就认识到地表以下流体流动中裂隙的作用，特别是近三十年来在油气藏工程、核废料的长期储存、边坡和坝肩稳定分析、地下洞室开挖稳定分析等领域裂隙的导体功能成为一个非常重要的问题。在水文地质领域研究地下水流动和污染物运移也会经常遇到断裂的地质结构，在这样的地质体里，大部分水流是通过裂隙传输的，因此，地下水流动和污染物运移由这些裂隙的特征所决定。裂隙的特征包括裂隙方向、长度，渗透性以及裂隙的密度和连通性。

离散裂隙网络模型(discretefracturenetwok，dfn)是研究裂隙岩体渗流最为有效的手段之一，dfn模型考虑的是不渗透或低渗透性岩石的裂隙岩体，而对于基质渗透性较大的裂隙岩体则在此基础上利用等效连续介质模型或双重介质模型来模拟。通常在裂隙介质中研究流体流动有三大主题：分别是连通性、流动性和扩散性，连通性排在第一，因此，裂隙的连通性是进行流体流动和溶质运移的前提和基础，裂隙的连通性对裂隙岩体的渗流有很大的影响。为了获得有效的渗透率样本值而避免无谓的计算采用了编程实现简单但效率相对较低的“染色”算法来研究多孔介质的连通，另外本行业对裂隙方向为均匀分布的等长裂隙和非等长裂隙的连通特性进行了研究。裂隙岩体渗流与应力耦合问题目前是岩石力学与岩石工程领域的热点问题之一，对于单一裂隙渗流–应力耦合模型目前试验研究成果较多，但对于多裂隙渗流–应力耦合模型，特别是包括随机裂隙网络建模及裂隙网络渗流路径识别等，研究较少，而这些成果都可以用于改进和扩展离散裂隙网络渗流–应力耦合模型，并应用于工程实践中。

技术实现要素：

针对上述现有技术存在的问题，本发明提供一种基于图论的裂隙网络连通性及渗流计算的方法，能模拟岩体多裂隙的渗流情况，并通过得出的渗流结果，为后续的工程实践提供理论支撑。

为了实现上述目的，本发明采用的技术方案是：该种基于图论的裂隙网络连通性及渗流计算的方法，具体步骤为：

a、二维裂隙网络裂隙迹长的确定：在二维裂隙网络中裂隙以迹线线段表示，裂隙迹长采用负指数分布或对数正态分布；

b、二维裂隙网络裂隙方向的确定；裂隙方向采用fisher分布，fisher分布的概率密度函数为：

式中：φ为与平均方向的偏差，k为fisher常数或分散因子，对每一条要产生的裂隙，随机地从fisher分布中取出的φ值(正值或负值)，加到平均裂隙方向θ即得到裂隙方向；

c、生成二维裂隙网络图：裂隙由中心点(x0,y0)、迹长l、θ角3个参数确定，迹线线段的端点坐标为

在设定的模拟区域内，迹线中心点服从均匀分布，中心点的个数由该区域裂隙的密度决定，裂隙条数依据其密度服从泊松分布，长度为迹长l，方向由θ角确定，θ角定义为自x轴逆时针旋转至迹线线段的角度；将生成的裂隙超出模拟区域边界的部分去除，通过montecarlo随机模拟方法生成该模拟区域的二维裂隙网络图；

d、生成二维裂隙网络的一级连通图：将二维裂隙网络图视作无向图，以无向图顶点表示裂隙，两个顶点间是否相连表示两条裂隙是否相交即是否有节点，通过无向图的邻接矩阵建立二维裂隙网络图的数据结构；将与模拟区域的左侧边界相交的裂隙组a找出，然后找寻出与上述裂隙组a相交的裂隙组b，重复上述过程，直至找出与模拟区域的右侧边界相交的裂隙组n，最终得出模拟区域的左侧边界、裂隙组a、裂隙组b…裂隙组n及模拟区域的右侧边界组成一级连通路径的各个裂隙，将无节点的裂隙删除，得出二维裂隙网络的一级连通图；

e、生成二维裂隙网络的二级连通图：将步骤d得出的一级连通图加入上下边界，计算裂隙之间以及裂隙与边界之间的交点，通过树裁剪算法依次删除裂隙网络中仅具有单个节点的裂隙，直至剩余裂隙至少具有两个节点，得出二级连通图；

f、计算二维裂隙网络的渗流情况：裂隙的导水系数采用对数正态分布，通过岩体裂隙渗流的立方定律和渗流连续性方程，建立裂隙网络的单相的、完全饱和的、稳定的渗流模型。左边界和右边界与裂隙相交的节点为边界节点，裂隙之间相交的节点为内部节点，边界节点的水头为已知，内部节点的水头为未知，采用裂隙网络的渗流模型在给定的边界条件下，即左右边界为定水头，上下边界为隔水边界，或左右边界为定水头，上下边界为从左边界到右边界的线性变化，完成内部节点水头的计算，得出内部节点的水头及整个模拟区域总的进出流量，进而根据节点的水头分布可绘制出水头等值线图，最终得出模拟区域的渗流场变化情况。

与现有技术相比，本发明通过二维裂隙网络中裂隙迹长的确定、二维裂隙网络中裂隙方向的确定、生成二维裂隙网络图、生成二维裂隙网络的一级连通图、生成二维裂隙网络的二级连通图及计算二维裂隙网络的渗流情况；可模拟出不同裂隙参数下，岩体裂隙的连通性及渗流特性，从而便于在工程实践中采集岩样得出该处岩体的裂隙参数分布，即可得出该处岩体的渗流情况，为工程实践的实施提供理论支撑。

附图说明

图1是本发明的整体流程图；

图2是本发明实施例中的二维裂隙网络图；

图3是本发明实施例中的一级连通图；

图4是本发明实施例中的二级连通图；

图5是本发明实施例中的水头等值线图。

具体实施方式

下面将对本发明作进一步说明。

实施例：如图所示，本发明包括具体步骤为：

a、二维裂隙网络中裂隙迹长的确定：在二维裂隙网络里裂隙以线段表示，裂隙迹长指岩体露头面和岩体结构面的交线长度，其能够代表结构面的规模并且在野外能够进行实际观测，裂隙的迹长概率分布一般服从负指数分布或对数正态分布；因此本发明在二维裂隙网络中裂隙以迹线线段表示，裂隙的迹长概率分布采用负指数分布或对数正态分布；

b、二维裂隙网络中裂隙方向的确定；岩体中的节理裂隙产状可能在一个或多个统计上占优的方向周围成组，裂隙系统往往由几类产状不同的节理组合而成，在大量测得节理产状后，同组裂隙在成因，分布特点，几何性质和工程性质等方面都具有相似性。本发明裂隙方向的概率分布采用fisher分布，fisher分布的概率密度函数为：

式中：φ为与平均方向的偏差，k为fisher常数或分散因子。高的k值意味着围绕平均方向更紧密的簇，该分布接近方差为1/k的正态分布，低的k值意味着围绕平均方向具有较宽的簇，当k＝0时，相当于均匀分布。对每一条要产生的裂隙，随机地从fisher分布中取出的φ值(正的或负的)，加到平均裂隙方向θ就得到了要产生的裂隙方向；裂隙中心点通常服从均均分布，裂隙条数依据其密度服从poisson随机过程。

c、生成二维裂隙网络图：裂隙由中心点(x0,y0)、迹长l、θ角3个参数确定，迹线线段的端点坐标为

在设定的模拟区域内，迹线中心点服从均匀分布，中心点的个数由该区域裂隙的密度决定，裂隙条数依据其密度服从泊松分布，长度为迹长l，方向由θ角确定，θ角定义为自x轴逆时针旋转至迹线线段的角度；将生成的裂隙超出模拟区域边界的部分去除，通过montecarlo随机模拟方法生成该模拟区域的二维裂隙网络图；

d、生成二维裂隙网络的一级连通图：将二维裂隙网络图视作无向图，以无向图顶点表示裂隙，两个顶点间是否相连表示两条裂隙是否相交即是否有节点，通过无向图的邻接矩阵建立二维裂隙网络图的数据结构；将与模拟区域的左侧边界相交的裂隙组a找出，然后找寻出与上述裂隙组a相交的裂隙组b，重复上述过程，直至找出与模拟区域的右侧边界相交的裂隙组n，最终得出模拟区域的左侧边界、裂隙组a、裂隙组b…裂隙组n及模拟区域的右侧边界组成一级连通路径的各个裂隙，将无节点(即无交点)的裂隙删除，得出二维裂隙网络的一级连通图；具体为：图2有24条裂隙(22条裂隙外加模拟区域的左右两条边界)相当于有24个顶点的无向图，根据平面内线段相交判断两裂隙是否相交，即无向图两顶点是否相连，建立二维裂隙网络图(见图2，表示为无向图g)的邻接矩阵：

a＝a(g)＝(aij)24×24(3)

其中，如果顶点vi邻接顶点vj，则aij＝1，否则aij＝0，于是24个顶点的无向图与对角线元素为零的24阶对称二元矩阵之间是一一对应的关系(见表1的a矩阵)；

表1二维裂隙网络图的数据结构(a矩阵)

找出a矩阵第一列为1的行数，9、13、14、23，共4行，也就是4条裂隙与左边界相交，与图2相对应，分别是裂隙9、13、14、23与左边界相交，首先从9行出发，寻找与9行相交的列为4列，再寻找与4列相交的行为3行，与3行相交的列有8、12、18、22、24，其中24列为右边界，说明本条路径连通左边界和右边界。继续寻找其它分支是否连到本条路径，退回到22列，与15行相交，说明裂隙22与15相交，图2上显示裂隙22与15相交，该分支再无其它路径，于是退回到18列，与18列相交的行有13、19、21，说明裂隙18与13、19、21相交，图2显示裂隙18与13、19、21相交，如此寻找下去，直至退回到4列，寻找出所有与9行相交的裂隙，本路径搜索结束，接下来用同样的方法对13行进行搜索，如此循环往复，直至把所有贯通左右边界的路径全部搜索出，从而提取出主干网(backbone)，该算法采用递归算法，称为一级连通图(见图3)。

e、生成二维裂隙网络的二级连通图：将步骤d得出的一级连通图加入上下边界(25线为下边界、26线为上边界)，计算裂隙之间以及裂隙与边界之间的交点，通过树裁剪算法依次删除裂隙网络中仅具有单个节点的裂隙，直至剩余裂隙至少具有两个节点，得出二级连通图；具体为：这些与主干网仅有1个节点的裂隙没有再连到别的裂隙或边界，为不导水的裂隙并进行删除，与主干网仅有1个节点的裂隙包括11，15，19，首先删除。在删除这些裂隙后，就会出现原来与主干网有2个节点的裂隙，现在变成1个节点，成为不导流裂隙如裂隙12，如此的树裁剪一直重复进行直到主干网上所有的裂隙至少有2个节点，意味着它们都是导流的，得到二级连通图(见图4)，用于裂隙渗流计算，表2记录下裁剪过程中裂隙与边界的连通性。

表2二维裂隙网络的连通矩阵

注：裂隙1是左边界，裂隙24是右边界，在一级连通图加下边界25线，上边界26线

f、计算二维裂隙网络的渗流情况：裂隙的导水系数采用对数正态分布，边界左侧与裂隙相交的节点为边界节点a，边界右侧与裂隙相交的节点为边界节点b，裂隙之间相交的节点为内部节点，边界节点a、b的水头为已知，内部节点的水头采用立方定律和渗流连续性方程建立的裂隙网络渗流模型进行求解，完成内部节点水头的计算，最终得出内部节点的水头及整个模拟区域总的进出流量。具体为：navier-stokes方程对层流、不可压缩流体在两个光滑平行板间流动的求解表达式为：

式中：q为单宽流量(m³/s/m)；μ为流体的动力粘滞系数(n·s/m²)；b为两个平行板间的距离或裂隙孔径(m)，为水力梯度；

则裂隙导水系数tf为：

于是

其中：c为水力传导率(m²/s)。

对于交叉裂隙水流运动，节点j质量守恒方程为：

式中，qij为节点i流向节点j的流量，ej为节点j的源汇项(如入渗量、开采量)。

如忽略源汇项，式(6)代入式(7)得：

式中，cij为节点i与j之间水力传导率，△hij＝hi-hj，其中hi、hj分别为节点i和节点j的水头(m)；

以两条相交裂隙水流为例，裂隙水流路径为具有五个相互连通的节点，假设基质不透水，裂隙中的流体是连续且不可压缩的，且连通节点处没有水头损失，那么在稳定流条件下的水量平衡方程为：

q15+q25+q35+q45＝0(9)

于是根据式(8)有

c15h1-c15h5+c25h2-c25h5+c35h3-c35h5+c45h4-c45h5＝0(10)

因此

于是，节点i和节点j水位关系为：

或

式中，其中，

假设模拟区域是在一个x-y的二维平面，则各点位置水头相同，则总水头就是压力水头；对每一个内部节点均可写出如公式(12)的一个方程。于是有n个内部节点就可列出n个方程，求解这个方程组，就可得到各个内部节点水头；

在选定的模拟区域内，其边界水头、水力梯度和裂隙导水系数确定后，可对裂隙网络进行渗流计算。本实施例设定左右边界为定水头边界，上下边界为隔水边界，水力梯度为0.2m/m，裂隙导水系数的均值为10^-6m²/s，标准差为10^0.2m²/s，基于公式(12)为每一个内部节点写出一个方程，如图4所示，图4包含17个节点，其中包含5个边界节点，12个内部节点，于是可写出12个方程，表3列出了12个方程组成的方程组的系数矩阵，比如对节点10，方程为：d8,10h8+d9,10h9-h10+d15,10h15+d17,10h17＝0。系数0表明该对节点间没有连通，如以表3中的下半部分内部水头的系数矩阵设为b矩阵，上半部分边界水头的系数矩阵设为b矩阵，边界水头矩阵为h0＝[h1h2h3h4h5]^t，则内部节点水头的求解就化为求解矩阵b'x＝-b'h0，利用matlab程序求解出x，即内部节点的水头(见表4)，根据节点的水头分布可绘制出水头等值线图(见图5)，流经整个模拟区域左边界的进水流量为qlb＝1.7706×10^-7m³/s/m，右边界的出水流量为qrb＝1.7706×10^-7m³/s/m，从而得出模拟区域的渗流场情况。

表3图4中12个内部节点流量平衡方程组成的系数矩阵

表4图4中的节点水头