[发明专利]基于裂隙网络连通性的裂隙流的MODFLOW模拟方法

权利要求书

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渗流计算和参数分析。