一种地下水大降深疏干下的全有效网格单元流场模拟方法与流程

本发明属于地下水系统数值仿真，尤其涉及一种地下水大降深疏干下的全有效网格单元流场模拟方法。

背景技术：

1、在过去数十年，包括人口增长、工业化加速发展、气候变化在内的多种因素导致全球地下水开采量明显增加，造成地下水含水层水位迅速下降(即“大降深”)甚至疏干(含水层中没有自由重力水面)，具有大降深疏干特点的地下水系统愈发常见。大降深疏干地下水系统引发了地面沉降、海水入侵、植被退化等一系列生态环境问题；准确模拟大降深疏干条件下的地下水流场对于科学合理地管理地下水系统的补给与排泄、恢复地下水位、修复相关生态环境问题具有重要支撑作用，已成为当前水资源保护和管理工作的重要支撑内容。

2、在地下水系统数值仿真中，往往将地下水含水层离散成具有指定行、列、层数的网格单元系统进行模拟，每个网格单元代表一定范围内的地下水含水层。目前，国际代表性地下水系统数值仿真模型如modflow对大降深疏干条件下的地下水系统的模拟能力不强。大降深疏干条件下地下水系统模拟的主要特点在于会出现大量被疏干的地下水系统数值仿真网格单元(以下简称“网格单元”，图1)。在经典的modflow-2005模型之中，疏干单元的出现会引发两方面问题，一是modflow-2005将疏干单元处理为无效网格单元，地下水的流动会在疏干单元处中断，地下水的流动在疏干单元(含水层)处不连续，这与实际情况不相符合；同时对于疏干单元，也不再模拟作用于该疏干单元(含水层)之上的源汇项，因此当大量网格单元疏干时会导致地下水系统水量平衡、地下水水头模拟结果失真。尽管modflow-2005可采用试算方法将符合预设湿润条件的疏干单元重新湿润使之再次成为有效单元，一定程度上可以减轻上述两方面问题带来的影响，但疏干单元的存在仍然在很大程度上限制了地下水系统数值仿真模型的应用面和准确性。为克服modflow-2005对疏干单元的处理方式所引起的弊端，在modflow的最新版本modflow-nwt中采用了特殊的“上风格式”算法来计算网格单元间的渗流量，据此可保证地下水的流动不在疏干单元处中断，疏干单元仍然可以被处理为有效网格单元，作用于疏干单元之上的源汇项模拟不会被取消，但由于“上风格式”算法本身精度较低，因此也从另一个方面损害了地下水系统数值仿真模型的准确性。

技术实现思路

1、为解决上述技术问题，本发明提出了一种地下水大降深疏干下的全有效网格单元流场模拟方法，提升了地下水系统数值仿真在大降深疏干条件下的物理意义和准确性，拓宽了地下水系统数值仿真的应用面。

2、为实现上述目的，本发明提供了一种地下水大降深疏干下的全有效网格单元流场模拟方法，包括以下步骤：

3、对地下水系统进行网格化剖分，获取地下水系统数值仿真网格单元；

4、基于达西定律结合调和平均渗透系数和加权平均过水断面厚度，获取沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数；

5、基于达西定律结合调和平均越流系数，获取沿垂向上地下水系统数值仿真所述第一网格单元和第三网格单元之间的水力传导系数；

6、基于水量平衡原则，获取地下水系统数值仿真所述第一网格单元的贮水项；

7、基于所述沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数、所述沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数及所述地下水系统数值仿真所述第一网格单元的贮水项，结合所述第一网格单元的水量平衡关系，构建大降深疏干条件下描述地下水运动的模型，迭代获取地下水系统各数值仿真所述第一网格单元的水头，通过图像化处理后获取大降深疏干条件下的流场。

8、可选的，获取所述地下水系统数值仿真网格单元的方法包括：

9、将地下水含水层空间上离散为具有一定行数、列数、层数的数值仿真网格单元系统，每个地下水系统数值仿真网格单元代表一定范围内的地下水含水层。

10、可选的，所述调和平均渗透系数和所述加权平均过水断面厚度的获取方法包括：

11、根据网格单元水平向的渗透系数、行间距、列间距、水头和顶底板高程，通过计算获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的所述调和平均渗透系数和所述加权平均过水断面厚度。

12、可选的，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的所述调和平均渗透系数的方法包括：

13、

14、其中，krave为沿水平方向上第一网格单元和第二网格单元之间的所述调和平均渗透系数；kri,j,k和kri,j+1,k分别为第一网格单元和第二网格单元沿行方向上的渗透系数；δrj+1/2为第一网格单元和第二网格单元中心点之间的距离；δrj和δrj+1分别为第j列、第j+1列网格单元的行间距。

15、可选的，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的加权平均过水断面厚度的方法包括：

16、在含水层没有顶板高程限制条件时，所述第一网格单元和所述第二网格单元之间的平均过水断面厚度通过网格单元各自的过水断面厚度线性加权计算如下，

17、

18、其中,hti,j,k和hti,j+1,k分别为第一网格单元和第二网格单元各自的过水断面厚度；hi,j,k和hi,j+1,k分别为这两个网格单元的水头；maxb为这两个网格单元中相对较高的底板高程；htave为第一网格单元和第二网格单元之间的加权平均过水断面厚度；

19、在含水层有顶板高程限制条件时，则网格单元间的加权平均过水断面厚度计算如下：

20、

21、其中，topi,j,k和topi,j+1,k分别为第一网格单元和第二网格单元各自的顶板高程；boti,j,k和boti,j+1,k分别为第一网格单元和第二网格单元各自的底板高程。

22、可选的，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的水力传导系数的方法包括：

23、

24、其中，cri,j+1/2,k为沿行方向上第一网格单元和第二网格单元之间的水力传导系数；δci为第i行网格单元的列间距；hup为第一网格单元和第二网格单元中水头相对较高的网格单元的水头。

25、可选的，获取沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数的方法及计算包括：

26、基于所述第一网格单元和所述第三网格单元之间的调和平均越流系数，及所述第一网格单元和所述第三网格单元的行间距、列间距，通过计算获取沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数；

27、cvi,j,k+1/2＝vci,j,k+1/2·δrj·δci

28、其中，cvi,j,k+1/2为沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数；vci,j,k+1/2为第一网格单元和第三网格单元之间的调和平均越流系数。

29、可选的，获取地下水系统数值仿真所述第一网格单元的贮水项的方法包括：

30、scii,j,k＝ui,j,k·δrj·δci

31、其中，scii,j,k为地下水系统数值仿真所述第一网格单元的贮水项，ui,j,k为贮水系数。

32、可选的，所述第一网格单元的水量平衡关系包括：

33、基于所述第一网格单元相邻的六个网格单元流入所述第一网格单元的水量与作用在所述第一网格单元之上的源汇项之和，去掉所述第一网格单元流出到所述六个网格单元的水量，获取所述第一网格单元的水量平衡关系。

34、可选的，构建大降深疏干条件下描述地下水运动的模型包括：

35、

36、其中，下标(i,j,k-1)，(i,j,k+1)，(i+1,j,k)，(i-1,j,k)，(i,j-1,k)，(i,j+1,k)分别为第一网格单元的上、下、前、后、左、右侧的相邻网格单元；cr、cc、cv分别为沿行、列、层方向上网格单元之间的水力传导系数；h为水头；上标n为本次迭代，上标n-1为上次迭代；pi,j,k为作用于第一网格单元之上的与水头有关的源汇项相关系数；qi,j,k为作用于第一网格单元之上的流量源汇项；scii,j,k为第一网格单元的贮水项；δt为模拟时段的时间长度；δh为模拟时段内的水头变化量。

37、本发明技术效果：本发明公开了一种地下水大降深疏干下的全有效网格单元流场模拟方法，首先对地下水系统进行网格化剖分，得到地下水系统数值仿真网格单元；然后根据达西定律，结合调和平均渗透系数、加权平均过水断面厚度、调和平均越流系数计算地下水系统数值仿真网格单元间的水平向和垂向水力传导系数，以表征含水层中地下水侧向、垂向流动的能力；然后根据水量平衡原理计算网格单元的贮水项；最后在上述计算结果的基础上根据水量平衡原理构建大降深疏干条件下的描述地下水运动的差分方程组，求解得到大降深疏干条件下的模拟流场。和现有技术相比，地下水的流动在疏干单元处不会被人为中断，疏干单元仍被处理为有效网格单元，作用于疏干单元之上的源汇项模拟不会被取消，同时网格单元间渗流量算法拥有很高的计算精度，因此提升了地下水系统数值仿真在大降深疏干条件下的物理意义和准确性，拓宽了地下水系统数值仿真的应用面。