CN109977353B

CN109977353B - 基于流网的非均质含水层优势流路径识别方法

Info

Publication number: CN109977353B
Application number: CN201910269272.4A
Authority: CN
Inventors: 南统超; 吴吉春; 王远坤
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2019-04-04
Filing date: 2019-04-04
Publication date: 2023-01-24
Anticipated expiration: 2039-04-04
Also published as: CN109977353A

Abstract

本发明公开了一种基于流网的非均质含水层优势流路径识别方法，该方法首先建立流管面积或体积与水流移动时间之间的定量关系；将已有的地下水流模拟程序用于非均质介质中流函数计算；通过计算流管面积或体积，高效、准确地识别时间最短的移动路径。与已有的粒子追踪法和图论技术相比，本发明提供的方法效率更高、准确性更好。

Description

基于流网的非均质含水层优势流路径识别方法

技术领域

本发明属于水力学技术领域，具体指代一种基于流网的非均质含水层优势流路径识别方法。

背景技术

目前，地下水污染问题得到越来越多的关注；地下水溶质运移难点之一在于地下水流场中优势流的存在，即：占含水层小部分体积、导水能力强的优势通道，常常可以传输大量的水流和溶质，甚至显著改变含水层行为。在自然界中普遍存在的含水层非均质性，是导致了地下水流场非均一性和优势流的主要因素。基于含水层非均质性识别地下水流场中的优势流路径，是地下水运移模拟与水质监控的重要部分。目前优势流路径识别与分析的主要手段是传统的粒子追踪法(PT)，以及近年新出现的基于图论的最小阻力路径法(MHR)。传统PT法在流场中释放大量粒子，通过追踪这些粒子在对流和弥散作用下的运动，分析水流运动不均匀性。该方法可靠性好，但主要问题在于计算量很大，需要大量时间及存储空间，而且根据粒子当前位置计算下一时刻位置容易引入不断积累的误差。MHR法以类似于导航寻路的算法寻找整体阻力最小的路径，计算效率较高，但其不考虑水动力条件和边界条件，容易导致不符合水动力规律的搜寻结果，因此该方法的可靠性还无法令人满意。

发明内容

针对于上述现有技术的不足，本发明的目的在于提供一种基于流网的非均质含水层优势流路径识别方法，通过利用已有的工具，便捷、高效且可靠地识别非均质含水层中的优势流路径，为有效刻画自然条件下污染物监控提供有力工具。

为达到上述目的，本发明采用的技术方案如下：

本发明的一种基于流网的非均质含水层优势流路径识别方法，包括步骤如下：

(1)根据所要研究的非均质含水层，确定渗透系数K的分布；

(2)对上述渗透系数K取倒数，得阻力场R分布，即R＝1/K；

(3)根据含水层条件，指定流函数Ψ的边界条件；

(4)计算含水层内流函数Ψ的分布；

(5)根据计算得到的流函数Ψ分布，采用等高线法计算一系列等差流线(即相邻流线的流函数之差相等)，流线包裹的空间即为流管；

(6)计算各流管的面积或体积；

(7)对各流管的面积或体积进行排序，面积或体积最小的流管对应最短的移动时间，即为优势流路径所在。

进一步地，所述步骤(3)中流函数Ψ与步骤(2)中阻力场R的关系满足如下方程：

其中，x为空间向量，该方程对应的边界条件需要指定如下：按照区域流场状况，上下游分别为流函数零坡度边界

其他平行于主流方向的边或面为恒定流函数边界Ψ(x)＝C，C为常数。

进一步地，所述步骤(4)中利用地下水流模拟软件(如MODFLOW、FEWATER等)在非均质含水层中求解上述流函数Ψ与阻力场R的关系方程，得到流函数的空间分布。

进一步地，所述步骤(5)利用等高线计算方法(如MATLAB中的contour函数等)计算若干等差流函数值的在空间上的分布，即等差流线。

进一步地，所述步骤(6)中使用polyarea函数计算流管的面积或体积。

进一步地，所述步骤(6)中计算的流管面积或体积，正比于流管内的渗流运动过程的时间t，即：

其中，ds为沿着流管方向的微小距离，v为达西流速的大小，w(s)为流管的宽度(三维情况下为截面积)，|Ψ_i+1-Ψ_i-1|为构成流管的流函数之差，A为该流管的总面积(三维情况下为总体积)。

本发明的有益效果：

1、本发明采用地下水模拟软件直接计算流函数，再通过流函数得到流线，不涉及大量计算，精度和效率均很高。

2、本发明既考虑了含水层非均质性，又能兼顾水动力条件，得出的优势流路径是流场中的一条流线，因此比MHR法的结果更符合物理规律、更可靠。

3、本发明各步骤涉及的细节技术可以用已有的工具解决，实现简单、操作便捷。

附图说明

图1为本发明的执行步骤示意图。

图2为流函数Ψ的边界条件示意图。

图3为实施例1中含有低渗透墙体的含水层示意图。

图4为实施例1中用本发明的SF方法得到的优势流路径结果示意图。

图5为实施例1中PT、MHR及SF方法得到的三种优势流路径结果与PT法中的粒子锋面及最先抵达处对比示意图。

图6为实施例2中含有低渗透墙体的含水层示意图。

图7为实施例2中PT、MHR及SF方法得到的优势流路径结果与PT法中的粒子锋面及最先抵达处对比示意图。

图8为实施例2中PT、MHR及SF方法得到的优势流路径结果与流场流线对比示意图。

图9为实施例3中对数渗透系数(1nK)空间分布示意图。

图10为实施例3中PT、MHR及SF方法得到的优势流路径结果与PT法中的粒子锋面及最先抵达处对比示意图。

图11为实施例3中PT、MHR及SF方法得到的优势流路径结果与流场流线对比示意图。

图12a为100个算例中观察到的第一种典型情况示意图。

图12b为100个算例中观察到的第二种典型情况示意图。

图12c为100个算例中观察到的第三种典型情况示意图。

图12d为100个算例中观察到的第四种典型情况示意图。

图13为100个算例中SF和MHR相对于PT路径的平均距离d的整体比较结果图。

图14为100个算例三种方法的耗时比较结果图。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

参照图1所示，本发明的一种基于流网的非均质含水层优势流路径识别方法，包括步骤如下：

(1)根据所要研究的非均质含水层，确定渗透系数K的分布；

(2)对上述渗透系数K取倒数，得阻力场R分布，即R＝1/K；

(3)根据含水层条件，指定流函数Ψ的边界条件；

(4)计算含水层内流函数Ψ的分布；

(6)计算各流管的面积或体积；

参照图2所示，所述步骤(3)中流函数Ψ与步骤(2)中阻力场R的关系满足如下方程：

以下对粒子追踪法简称PT，基于图论的最小水力阻力法简称MHR，本发明的基于流函数的算法简称SF。

实施例1：含有两个低渗透墙体的均质含水层

研究对象为含有两个低渗透墙体的均质承压含水层，区域平均水流由左到右，主体与低渗墙体的渗透系数比为K/K_b＝10³(如图3所示)，孔隙度

为常数。采用本发明的SF方法，计算等差流线和各流管面积A(如图4所示，图中显示了等差流线的空间分布，即相邻流线的流函数之差相等；右侧展示了流管的A值，最小A值以空心圆标示，接近含水层对称中轴)，最小的A对应的流管即为优势流路径。所得的SF路径为一平直通过墙体空隙的流线，与常识相符。为了比较本发明与PT、MHR方法并进行进一步分析，图5将PT、MHR及SF方法得到的三种优势流路径结果与PT法中的粒子锋面及最先抵达处进行了对比。结果表明：SF路径(图5中虚线)与PT路径(图5中实线)几乎完全一致，均指向粒子最先抵达处；MHR路径(图5中点线)也为一通过墙体的平直线。实际上所有通过空隙的平直线对MHR算法来说都是优势路径，然而真正的最快路径只有PT路径那一条(当然，算得的PT路径本身也有一定的误差)。但结合水动力条件可知，除了对称轴y＝100和上下边界以外的流线都是弯曲的，MHR方法无法识别这些路径之间的差别，也不考虑这些路径是否在水动力学上是否真实可行的。尽管本例中三种方法的最终结果差异不太大，但在复杂含水层中SF与MHR的差异会变得更加显著。

实施例2：含有多个低渗透墙体的均质含水层；

与前例相似，考虑一个均值承压含水层中低渗透墙体的情况，不同的是低渗墙体更多、分布更复杂(如图6所示)。同样，用PT、MHR、SF三种方法寻找优势流路径，并与PT方法中的粒子锋面及最先抵达处比较(见图7)。结果表明，本发明的SF路径(虚线)与PT路径(实线)几乎完全一致；MHR路径则有显著差异，体现为规避低渗墙体的最短路径。图8进一步展示了本例中水流场的流线分布以及各流管的A值，图8右侧展示了各流管的A值(最小A值以空心圆标示)。表明PT路径和SF路径与水流状态相符，而MHR路径切割多条流线，完全不考虑具体水流因素。这进一步表明本发明的SF路径的合理性与可靠性。

实施例3：多高斯对数渗透系数场中的优势流

考虑二维矩形承压含水层，对数渗透系数lnK满足多高斯分布假设，lnK的均值为0，方差为1.5，相关长度为50m。用序贯高斯模拟生成一个随机lnK样本(如图9所示)。同样考虑从左向右的区域水流，三种方法寻找的优势流路径与PT方法中粒子锋面、最先抵达处的对比见图10。结果表明，SF路径与PT路径基本一致，均指向粒子最先抵达处；MHR路径则剧烈摇摆，与PT路径有显著差别，穿过多个粒子锋面，右边终点也与最先抵达路径不一致。说明本例中SF路径远比MHR路径更接近真正的优势流路径。图11比较了三种方法获得的路径与流线分布的关系，右侧展示了各流管的A值(最小A值以空心圆标示)，结果表明PT路径和SF路径与流场相符，而MHR路径反复斜跨多条流线，不符合水动力规律。因此，SF路径更准确，也更合理。

实施例4：100个算例的整体比较

由于含水层具体的非均质分布可能影响各方法的效果和时间消耗，为了避免因个别的渗透系数分布影响对方法的判断，用前述的序贯高斯模拟随机生成了100个对数渗透系数lnK场，对每个lnK场都用三种方法进行优势流路径识别，对整体效果做比较。在这100个平行实验中，观察到四种典型结果，即：(a)含水层存在唯一显著优势通道，三种方法结果几乎完全一致(图12a)；(b)含水层存在唯一显著优势通道，SF路径与PT路径基本吻合，MHR路径偏差很大(图12b)；(c)含水层存在多个优势通道，SF路径与其他两种不同但与某一通道吻合(图12c)；(d)含水层存在多个优势通道，SF路径与PT路径吻合很好，MHR路径显著不同且不与任一通道吻合良好(图12d)。

为了定量比较SF路径或MHR路径与作为参照的PT路径的差异，引入平均路径距离d，计算方法如下：

其中，r(i)≡min_{j＝1，...，N′}|X_tes(i)-x_PT(j)|，X_tes(i)，i＝1，...，N为均匀离散的待测路径(即MHR路径或SF路径)，X_PT(j)，j＝1，...，N′为作为参照的PT路径。r(i)表示待测路径的一点到PT路径X_PT的最小距离。100个随机算例中的平均路径距离d如图13所示。结果表明SF路径所得d的均值、中位数、最小值、方差、四分位数等都比MHR路径d值的相应统计量更小。但SF路径d值的最大值比MHR路径d值最大值更大。经粒子锋面位置验证，这是由于流场中存在多个优势流通道且SF路径反映的是其他优势流路径。多个优势流通道的存在是SF路径某些d值较大的原因。而MHR路径常常不与任一优势流通道吻合反而常常斜跨多个可能的通道，导致其d值偏大但又比走向基本平行的SF路径d值稍小。

图14是100个算例中三种方法的耗时比较。显然，SF方法耗时最少，传统PT方法耗时最多，且远多于其他两种方法。

从这些算例中可以看出，本发明提供的SF方法获得的优势路径，基本与传统的PT方法提供的结果一致，存在多个优势流路径时，于其中一个路径吻合较好，考虑了水动力条件，比MHR方法的结果更合理。同时，SF方法的计算效率远高于PT方法，也高于MHR方法。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims

1.一种基于流网的非均质含水层优势流路径识别方法，其特征在于，包括步骤如下：

(1)根据所要研究的非均质含水层，确定渗透系数K的分布；

(2)对上述渗透系数K取倒数，获得阻力场R分布，即R＝1/K；

(3)根据含水层条件，指定流函数Ψ的边界条件；

(4)计算含水层内流函数Ψ的分布；

(5)根据计算得到的流函数Ψ分布，采用等高线法计算一系列等差流线，流线包裹的空间即为流管；

(6)计算各流管的面积或体积；

(7)对各流管的面积或体积进行排序，面积或体积最小的流管对应最短的移动时间，即为优势流路径所在；

所述步骤(3)中流函数Ψ与步骤(2)中阻力场R的关系满足如下方程：

其他平行于主流方向的边或面为恒定流函数边界Ψ(x)＝C，C为常数；

所述步骤(4)中利用地下水流模拟软件在非均质含水层中求解上述流函数Ψ与阻力场R的关系方程，得到流函数的空间分布；

所述步骤(6)中使用polyarea函数计算流管的面积或体积；

所述步骤(6)中计算的流管面积或体积，正比于流管内的渗流运动过程的时间t，即：

其中，ds为沿着流管方向的微小距离，v为达西流速的大小，w(s)为流管的宽度，|Ψ_i+1-Ψ_i-1|为构成流管的流函数之差，A为该流管的总面积或总体积。

2.根据权利要求1所述的基于流网的非均质含水层优势流路径识别方法，其特征在于，所述步骤(5)利用等高线计算方法计算若干等差流函数值的在空间上的分布，即等差流线。