CN108399290B - 基于裂隙网络连通性的裂隙流的modflow模拟方法 - Google Patents

Abstract

本发明公开了一种基于裂隙网络连通性的裂隙流的MODFLOW模拟分析方法，该方法通过利用Monte Carlo随机模拟技术生成裂隙网络，并且建立DFN模型，然后通过DFN和SC方法进行渗流分析。本发明在连续—非连续介质模型的构架内，将单条裂隙刻画到有限差分网格，不但有效地保存下了裂隙流，同时还能应用传统的地下水流模拟软件如MODFLOW进行渗流计算，为裂隙岩体的渗流计算提供了一种简单、实用的计算方法。

Description

基于裂隙网络连通性的裂隙流的MODFLOW模拟方法

技术领域

本发明属于裂隙岩体的渗流特性及其数值模拟，具体涉及一种基于裂隙网络连通性的裂隙流的MODFLOW模拟方法。

背景技术

20世纪70年代以来，裂隙岩体的渗流特性及其数值模拟一直是国内外学者研究的热点问题，随着计算机技术的发展，数值方法越来越多地用于岩石工程问题。省时、省力、灵活、经济是数值计算相比于理论分析和试验手段所具有的优势，而有效模拟岩体的结构特性，却是岩石工程数值计算的前提和挑战。天然岩体中存在大量的孔隙或裂隙，这些孔隙或裂隙不但改变了岩体的力学性质，也严重影响着岩体的渗透特性；它们不但具有非连续性、非均质性、各向异性，还具有多尺度效应。

等效连续介质模型(Stochastic Continuum,SC)假定裂隙岩体为服从达西流的均质多孔介质，不考虑单个裂隙的几何结构；而离散裂隙网络模型(Discrete FractureNetwork,DFN)中，考虑的是具有低渗透性基质的裂隙岩体，首先对研究区域岩体裂隙的现场测量和统计的基础上生成随机裂隙网络，然后运用渗流连续性方程和立方定律对裂隙岩体进行渗流求解；连续—非连续介质模型一般将岩体视为等效连续介质(不考虑裂隙网络)和裂隙网络的叠加，分别对岩体渗流进行求解，然后再进行叠加。等效连续介质模型的优点是模型简单、理论成熟、容易求解，适合于研究大范围、裂隙分布密集的岩体渗流，而对于岩体典型单元体很大甚至不存在的研究区域或者需要研究某个特定小范围、特定路径的渗流行为时，该模型却不适用，而离散裂隙网络模型考虑的是具有低渗透性基质的裂隙岩体，对于基质渗透性较大的岩体模拟是不适用的，而且研究尺度较大，往往导致求解困难；连续—非连续介质模型从原理上能够比较真实地刻画渗流形态，但裂隙网络和等效介质之间的水量交换往往难以确定。

发明内容

发明目的：针对上述现有技术存在的不足，本发明提供一种基于裂隙网络连通性的裂隙流的MODFLOW模拟方法。

技术方案：基于裂隙网络连通性的裂隙流的MODFLOW模拟方法，包括如下步骤：

(1)建立二维裂隙网络：通过Monte Carlo随机模拟技术和Matlab程序生成二维裂隙网络，包括生成裂隙网络原图、一级裁剪连通图、二级裁剪连通图和三级裁剪连通图；

(2)DFN渗流模拟：利用立方定律求解Navier-Stokes方程计算渗流，包括对于完全饱和的、单相的、不可压缩的层流，其计算表达式为：

其中，ΔH为裂隙端点的水头差；

(3)SC渗流模拟：将长度为l、导水系数为T_f的裂隙画到一个x-y的二维平面，使用地下水流数值模拟软件如MODFLOW计算渗流；

(4)DFN和SC渗流分析：以DFN渗流作为度量标准，分析不同网格分辩率对SC渗流的影响。

进一步的，所述的步骤(1)具体包括如下步骤：

(11)确定迹线参数，所述的参数包括裂隙中心点(x₀,y₀)、迹线长度l和迹线方向θ(自x轴逆时针旋转至迹线的角度)，建立迹线端点坐标(x，y)，所述的迹线端点坐标为：

端点坐标利用Matlab程序表示为：

x(1,:)＝p(1,:)+p(3,:)/2.*cos(p(4,:))

x(2,:)＝p(1,:)-p(3,:)/2.*cos(p(4,:))

y(1,:)＝p(2,:)+p(3,:)/2.*sin(p(4,:))

y(2,:)＝p(2,:)-p(3,:)/2.*sin(p(4,:))

其中：p(1,:)，p(2,:)表示迹线中心点的坐标，p(3,:)表示迹线长，p(4,:)表示迹线方向；

(12)生成二维裂隙网络的原图和一级裁剪连通图：对裂隙的迹线参数应用MonteCarlo随机模拟技术和Matlab程序生成二维裂隙网络原图，然后通过无向图的邻接矩阵建立二维裂隙网络的数据结构，利用递归算法搜索出所有贯通左右边界的路径，提取出主干网为一级裁剪连通图；

(13)生成二维裂隙网络的二级裁剪连通图：利用树裁剪算法删除一级裁剪连通图中只有单个交点的裂隙，重复进行直到主干网上所有的裂隙至少有2个交点，得到二级裁剪连通图；

(14)生成二维裂隙网络的三级裁剪连通图：二级裁剪连通图用于DFN渗流计算，但是对于SC网格渗流计算时，二级裁剪连通图的裂隙死端(fracturedead-end)会导致网格间错误的连通，应当删除，从而得到三级裁剪连通图。

进一步的，所述的步骤(2)具体步骤如下：

(21)对于完全饱和的、单相的、不可压缩的层流，通过立方定律求解Navier—Stokes方程，其计算表达式为：

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

为水力梯度；

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

则

式中，C_ij为节点i与j之间水力传导率，ΔH_ij＝H_i-H_j，其中H_i、H_j分别为节点i和节点j的水头(m)。

进一步的，所述的步骤(3)包括如下步骤：

(31)走向为x方向或y方向的裂隙刻画到单元格大小为Δ的有限差分网格上，该裂隙穿过每一个单元格的x方向或y方向的渗透系数估计为T_f/Δ；

(32)基于裂隙与网格单元格的对齐方式将裂隙分为两种类型，分别为类型I和类型II，所述的类型I裂隙为经过单元格两个相对的面进出单元格，所述的类型II裂隙指经过单元格并不相对的两个面进出单元格；

(33)如果在某个单元格仅有类型II的一条裂隙，其渗透系数K_MODFLOW分配到水平和垂直两个方向，如果在同一单元格里有多条类型II的裂隙，则类型II裂隙的K_MODFLOW平分到水平和垂直两个方向，然后每个方向的渗透系数相加，因此为每个单元格可计算出不同的各向异性比率。

更进一步的，所述的步骤(4)以DFN渗流作为度量标准，分析不同网格分辩率对SC渗流的影响。包括利用二维裂隙网络的三级裁剪连通图进一步细化裂隙网格，进行不同网格分辨率的SC渗流分析。

有益效果：本发明相比现有技术其显著的效果在于，本方法可以视为在连续—非连续介质模型的构架内，将单条裂隙刻画到有限差分网格，不但有效地保存下了裂隙流，同时还能应用传统的地下水流模拟软件如MODFLOW进行渗流计算，为裂隙岩体的渗流计算提供了一种简单、实用的计算方法。

附图说明

图1本发明的处理步骤及其流程图；

图2为本发明岩体中迹线的控制参数示意图；

图3为本发明实施例1裂隙网络的一次随机实现的二维裂隙网络原图；

图4为本发明四个顶点的无向图的示意图及其邻接矩阵；

图5为本发明实施例1的一级裁剪连通图；

图6为本发明实施例1的二级裁剪连通图；

图7为本发明实施例1的三级裁剪连通图；

图8为本发明五个连通节点形成的一个典型网络示意图；

图9为本发明沿网格的流动路径示意图(水平线表示类型I的裂隙网格，倾斜线表示类型II的裂隙网格)；

图10为本发明实施例2裂隙网络的一次随机实现的二维裂隙网络原图；

图11为本发明实施例2中的二级裁剪连通图；

图12为本发明实施例2中的三级裁剪连通图。

具体实施方式

为了详细的说明本发明公开的技术方案，下面结合具体实施例和说明书附图做进一步的阐述。

实施例1

本发明所提供的一种基于裂隙网络连通性的裂隙流的MODFLOW模拟方法，参见图1所示，主要包括以下步骤：

1、建立二维裂隙网络：通过Monte Carlo随机模拟技术和Matlab程序生成二维裂隙网络，包括生成裂隙网络原图、一级裁剪连通图、二级裁剪连通图和三级裁剪连通图。

在一个二维裂隙网络中，迹线通常表示为线段，由裂隙中心点(x₀,y₀)、迹线长度l、裂隙方向θ角3个参数确定，θ角定义为自x轴逆时针旋转至迹线的角度，如图2所示。

迹线端点坐标(x,y)为：

给定产生的裂隙数(如实施例1)或给出每一组裂隙的密度以确定该组产生的裂隙数(如实施例2)，裂隙中心通常服从均匀分布，裂隙长度和方向分别从各自的分布中随机选取，以p(1,:),p(2,:)表示迹线中心点的坐标，p(3,:)表示迹线长，p(4,:)表示迹线方向，利用Matlab程序裂隙的端点坐标表示如下：

x(1,:)＝p(1,:)+p(3,:)/2.*cos(p(4,:))；

x(2,:)＝p(1,:)-p(3,:)/2.*cos(p(4,:))；

y(1,:)＝p(2,:)+p(3,:)/2.*sin(p(4,:))；

y(2,:)＝p(2,:)-p(3,:)/2.*sin(p(4,:))；

裂隙与模型边界相交，将超出模型边界的裂隙裁减掉，利用Monte Carlo随机模拟技术生成二维裂隙网络原图，如图3所示，实施例1的裂隙统计参数见表1。

表1实施例1的裂隙统计参数

由于裂隙网络与图论的无向图具有相拟的特点，以无向图顶点表示裂隙，两个顶点间是否相连表示两条裂隙是否相交，于是通过无向图的邻接矩阵建立裂隙网络的数据结构，无向图的示意图如图4所示；根据图3建立的数据结构利用递归算法把图3中所有贯通左右边界的路径全部搜索出，从而提取出主干网，称为一级裁剪连通图，如图5所示；再利用树裁剪算法删除一级裁剪连通图中只有单个交点的裂隙，这些与主干网仅有1个交点的裂隙没有再连到别的裂隙或边界，认为是不导水的，需要删除，在删除这些裂隙后，就会出现有一些别的裂隙原来与主干网有2个交点，现在变成1个交点，成为不导流裂隙，同样也要删除，如此的树裁剪一直重复进行直到主干网上所有的裂隙至少有2个交点，意味着它们都是导流的，得到二级裁剪连通图见图6。

由于DFN渗流是基于连接节点处每条裂隙的路径，因此二级裁剪连通图可用于DFN渗流的计算，可是对于SC的网格渗流计算比如用MODFLOW,图6的裂隙死端会导致网格间错误的连通，应当删除，从而得到三级裁剪连通图如图7，用于网格渗流计算。

2、DFN渗流模拟

对于完全饱和的、单相的、不可压缩的层流，Navier-Stokes方程使用立方定律求解其表达式为：

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

为水力梯度。

则裂隙导水系数T_f为：

于是

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

对于交叉裂隙水流运动，如忽略源汇项，则节点j质量守恒方程为：

即

式中，C_ij为节点i与j之间水力传导率，ΔH_ij＝H_i-H_j，其中H_i、H_j分别为节点i和节点j的水头(m)。

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

Q₁₅+Q₂₅+Q₃₅+Q₄₅＝0 (7)

于是根据式(6)有

C₁₅H₁-C₁₅H₅+C₂₅H₂-C₂₅H₅+C₃₅H₃-C₃₅H₅+C₄₅H₄-C₄₅H₅＝0 (8)

因此

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

或

式中，

其中，

对每一个内部节点都可写出如公式(10)的一个方程。于是在实施例1有12个内部节点就可列出12个方程，利用Matlab程序求解该方程组，就可得到内部节点水头。

3、SC渗流模拟

SC渗流模拟是将裂隙网络随机生成的单个裂隙映射到连续渗流场中的有限差分网格上。我们将长度为l、导水系数为T_f的裂隙画到一个x-y的二维平面，流体在裂隙缝隙的单宽流量为：

其中ΔH为裂隙端点的水头差。

对于走向为x方向的裂隙画到单元格大小为Δ的有限差分网格上，该裂隙穿过每一个单元格的x方向的渗透系数估计为T_f/Δ；同理，对于走向为y方向的裂隙画到单元格大小为Δ的有限差分网格上，该裂隙穿过每一个单元格的y方向的渗透系数估计为T_f/Δ，如图9所示。

因为与x方向成一定角度的裂隙不再与网格对齐(见图9)，较长的流动路径(L sin(θ)+L cos(θ))包含了水平和垂直分量，导致水力梯度降低，因此单元格的渗透系数需要增大来补偿水力梯度以便保存下相同的流量。使用地下水流数值模拟软件如MODFLOW，渗透系数K表示为：

其中Δ为单元格大小。

基于裂隙与网格单元格的对齐方式将裂隙分为两种类型：类型I和类型II。类型I裂隙指经过单元格两个相对的面进出单元格，比如水平或垂直裂隙；而类型II裂隙指经过单元格并不相对的两个面进出单元格，比如倾斜裂隙。

如果在某个单元格仅有类型II的一条裂隙，其渗透系数K_MODFLOW分配到水平和垂直两个方向，如果在同一单元格里有多条类型II的裂隙，则类型II裂隙的K_MODFLOW平分到水平和垂直两个方向，然后每个方向的渗透系数相加，因此为每个单元格可计算出不同的各向异性比率。

裂隙的有限差分映射和渗透系数的修正在小网格分辨率取得了很好的流量估计(同DFN渗流比较)。但随着单元格分辨率的增大，差异会迅速增大。

4、DFN和SC渗流分析

对图6(DFN渗流)和图7(SC渗流)的边界条件设定为上下为隔水边界，左右为定水头边界(分别设为10m和5m)。其裂隙流为：Q_DFN＝2.9001×10^-7m³/sec/m，SC分辨率为0.05和0.1m的网格流分别为Q_GRID＝2.9394×10^-7m³/sec/m和Q_GRID＝3.1376×10^-7m³/sec/m，同DFN流的误差分别为1.36％和8.19％。表2列出了节点坐标和相应的DFN水头和SC两个不同分辨率的水头。SC细网格(如0.05m)内部节点水头和相应的DFN水头误差为0.00～0.01m，而SC较粗网格(如0.1m)内部节点水头和相应的DFN水头误差为0.01～0.13m。结果表明细网格的MODFLOW模拟获得了较好的结果，但是粗网格极大地高估了流量和水头。

表2实施例1裂隙节点坐标和DFN及SC模拟的水头

在图6中裁剪掉裂隙死端得到的用于SC渗流计算的图7中，粗网格可能导致缩短或延长流动路径，但细网格的流动路径总是趋近于真实的裂隙长度见表3。

表3实施例1的DFN和SC粗网格和细网格估计的流动路径长度(m)

实施例2

为了验证该结论，再次设计了较为复杂的裂隙网络，参见图10所示，实施例2的裂隙统计参数见表4，图11、图12分别为实施例2对应的二级、三级裁剪连通图。实施例2的边界条件同实施例1，其裂隙流Q_DFN＝1.6865×10^-6m³/sec/m，SC分辨率为0.05和0.1m的网格流分别为Q_GRID＝1.7269×10^-6m³/sec/m和Q_GRID＝1.8078×10^-6m³/sec/m，同DFN流的误差分别为2.40％和7.19％，SC细网格(如0.05m)内部节点水头和相应的DFN水头误差为0.00～0.08m，而SC粗网格(如0.1m)内部节点水头和相应的DFN水头误差为0.00～0.16m。作为比较，将网格的分辨率进一步缩小到0.02m，其Q_GRID＝1.710×10^-6m³/sec/m,与DFN流的误差仅为1.40％，MODFLOW模拟的水头非常接近DFN模拟，它们之间的误差为0.00～0.03m(见表5)。

表4实施例2的裂隙统计参数

表5实施例2裂隙节点坐标和DFN及SC模拟的水头

MODFLOW模拟在实施例1的大部分节点高估了水头，而在实施例2却是大部分节点低估了水头(表5中H(MOD,0.02),H(MOD,0.05)和H(MOD,0.1)分别减去H(DFN)得到的负值)。通过网格细化低估值极大地降低了。大的低估值发生在裂隙网络的裂隙分布稠密的节点处(如图12节点

和

)。

低估的水头表明由于粗网格拉长了裂隙长度从而导致邻近节点间的水力梯度低估，从表6可以看出这些拉长的裂隙长度。比如在粗网格中节点

被拉长0.08m，节点

是0.03m，节点

是0.06m，这些被拉长的裂隙均发生在裂隙分布稠密的节点处。其高估长度的累积效应大大低估了流量和水头。

表6实施例2的DFN和SC粗网格和细网格估计的流动路径长度(m)

与此相反的是邻近边界的节点

和

也处于稠密的裂隙分布区域，这些节点间的长度也是高估了，如节点

高估了0.04m，节点

和

均高估了0.05m，但是，这些节点间MODFLOW模拟的水头误差比接近研究区域中心稠密的区域却小得多，表明给定的边界条件约束了它们邻近节点的水头计算从而降低了它们的计算误差。

Claims (4)

1.基于裂隙网络连通性的裂隙流的MODFLOW模拟方法，其特征在于：包括如下步骤：

(1)建立二维裂隙网络：通过Monte Carlo随机模拟技术和Matlab程序生成二维裂隙网络，包括裂隙网络原图、一级裁剪连通图、二级裁剪连通图和三级裁剪连通图；

(2)DFN渗流模拟：利用立方定律求解Navier-Stokes方程并计算渗流，包括对完全饱和的、单相的、不可压缩的层流的渗流参数求解，其计算表达式如下：

其中，ΔH为裂隙端点的水头差；

(3)SC渗流模拟：将长度为l、导水系数为T_f的裂隙画到一个x-y的二维平面，使用地下水流数值模拟软件MODFLOW计算渗流；

(4)DFN和SC渗流分析：以DFN渗流作为度量标准，分析不同网格分辨率对SC渗流的影响；

所述步骤(3)具体步骤如下：

(31)将走向为x方向或y方向的裂隙画到单元格大小为Δ的有限差分网格上，该裂隙穿过每一个单元格的x方向或y方向的渗透系数估计为T_f/Δ；

(32)将基于裂隙与网格单元格的对齐方式将裂隙分为两种类型，分别为类型I和类型II，所述的类型I裂隙为经过单元格两个相对的面进出单元格，所述的类型II裂隙指经过单元格并不相对的两个面进出单元格；

(33)如果在某个单元格仅有类型II的一条裂隙，其渗透系数K_MODFLOW分配到水平和垂直两个方向，如果在同一单元格里有多条类型II的裂隙，则类型II裂隙的K_MODFLOW平分到水平和垂直两个方向，然后每个方向的渗透系数相加，因此为每个单元格计算出不同的各向异性比率。

2.根据权利要求1所述的基于裂隙网络连通性的裂隙流的MODFLOW模拟方法，其特征在于：所述的步骤(1)具体步骤如下：

(11)确定迹线参数：所述的参数包括裂隙中心点(x₀，y₀)、迹线长度l和迹线方向θ，建立迹线端点坐标(x，y)，所述的迹线端点坐标为：

端点在Matlab程序中坐标如下：

x(1,:)＝p(1,:)+p(3,:)/2.*cos(p(4,:))

x(2,:)＝p(1,:)-p(3,:)/2.*cos(p(4,:))

y(1,:)＝p(2,:)+p(3,:)/2.*sin(p(4,:))

y(2,:)＝p(2,:)-p(3,:)/2.*sin(p(4,:))

其中，p(1,:)、p(2,:)表示迹线中心点的坐标，p(3,:)表示迹线长，p(4,:)表示迹线方向；

(12)生成二维裂隙网络的原图和一级裁剪连通图：对裂隙的迹线参数应用MonteCarlo随机模拟技术和Matlab程序生成二维裂隙网络原图，然后通过无向图的邻接矩阵建立二维裂隙网络的数据结构，利用递归算法搜索出所有贯通左右边界的路径，提取出主干网为一级裁剪连通图；

(13)生成二维裂隙网络的二级裁剪连通图：利用树裁剪算法删除一级裁剪连通图中只有单个交点的裂隙，重复进行直到主干网上所有的裂隙至少有2个交点，得到二级裁剪连通图；

(14)生成二维裂隙网络的三级裁剪连通图：删除二级裁剪连通图的裂隙死端，得到三级裁剪连通图。

3.根据权利要求1所述的基于裂隙网络连通性的裂隙流的MODFLOW模拟方法，其特征在于：步骤(2)具体步骤如下：

(21)对于完全饱和的、单相的、不可压缩的层流，通过立方定律求解Navier-Stokes方程，其计算表达式如下：

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

为水力梯度；

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

则

式中，C_ij为节点i与j之间水力传导率，ΔH_ij＝H_i-H_j，其中H_i、H_j分别为节点i和节点j的水头(m)。

4.根据权利要求1所述的基于裂隙网络连通性的裂隙流的MODFLOW模拟方法，其特征在于：所述的步骤(4)以DFN渗流作为度量标准，包括利用二维裂隙网络的三级裁剪连通图进一步细化裂隙网格，进行不同网格分辨率的SC渗流计算和参数分析。

Images