CN116992783A - 一种地下水大降深疏干下的全有效网格单元流场模拟方法 - Google Patents
一种地下水大降深疏干下的全有效网格单元流场模拟方法 Download PDFInfo
- Publication number
- CN116992783A CN116992783A CN202310582904.9A CN202310582904A CN116992783A CN 116992783 A CN116992783 A CN 116992783A CN 202310582904 A CN202310582904 A CN 202310582904A CN 116992783 A CN116992783 A CN 116992783A
- Authority
- CN
- China
- Prior art keywords
- grid
- grid cell
- water
- groundwater
- numerical simulation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 187
- 238000004088 simulation Methods 0.000 title claims abstract description 91
- 238000000034 method Methods 0.000 title claims abstract description 42
- 239000003673 groundwater Substances 0.000 claims abstract description 73
- 230000035699 permeability Effects 0.000 claims abstract description 30
- 238000003860 storage Methods 0.000 claims abstract description 29
- 230000033001 locomotion Effects 0.000 claims abstract description 15
- 238000003384 imaging method Methods 0.000 claims abstract description 7
- 238000004364 calculation method Methods 0.000 claims description 22
- CIWBSHSKHKDKBQ-JLAZNSOCSA-N Ascorbic acid Chemical compound OC[C@H](O)[C@H]1OC(=O)C(O)=C1O CIWBSHSKHKDKBQ-JLAZNSOCSA-N 0.000 claims description 16
- ZZZCUOFIHGPKAK-UHFFFAOYSA-N D-erythro-ascorbic acid Natural products OCC1OC(=O)C(O)=C1O ZZZCUOFIHGPKAK-UHFFFAOYSA-N 0.000 claims description 8
- 229930003268 Vitamin C Natural products 0.000 claims description 8
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 claims description 8
- 235000019154 vitamin C Nutrition 0.000 claims description 8
- 239000011718 vitamin C Substances 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 239000006185 dispersion Substances 0.000 description 3
- 238000006467 substitution reaction Methods 0.000 description 3
- 238000009736 wetting Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 230000008093 supporting effect Effects 0.000 description 2
- 230000015556 catabolic process Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000001771 impaired effect Effects 0.000 description 1
- 230000009545 invasion Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- JTJMJGYZQZDUJJ-UHFFFAOYSA-N phencyclidine Chemical compound C1CCCCN1C1(C=2C=CC=CC=2)CCCCC1 JTJMJGYZQZDUJJ-UHFFFAOYSA-N 0.000 description 1
- 239000013535 sea water Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Computing Systems (AREA)
- Software Systems (AREA)
- Algebra (AREA)
- Computer Graphics (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种地下水大降深疏干下的全有效网格单元流场模拟方法，包括以下步骤：对地下水系统进行网格化剖分，获取地下水系统数值仿真网格单元；基于达西定律结合调和平均渗透系数和加权平均过水断面厚度，获取沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数；基于达西定律结合调和平均越流系数，获取沿垂向上地下水系统数值仿真所述第一网格单元和第三网格单元之间的水力传导系数；基于水量平衡原则，获取地下水系统数值仿真所述第一网格单元的贮水项；构建大降深疏干条件下描述地下水运动的模型，迭代获取地下水系统各数值仿真所述第一网格单元的水头，通过图像化处理后获取大降深疏干条件下的流场。
Description
技术领域
本发明属于地下水系统数值仿真技术领域，尤其涉及一种地下水大降深疏干下的全有效网格单元流场模拟方法。
背景技术
在过去数十年，包括人口增长、工业化加速发展、气候变化在内的多种因素导致全球地下水开采量明显增加，造成地下水含水层水位迅速下降(即“大降深”)甚至疏干(含水层中没有自由重力水面)，具有大降深疏干特点的地下水系统愈发常见。大降深疏干地下水系统引发了地面沉降、海水入侵、植被退化等一系列生态环境问题；准确模拟大降深疏干条件下的地下水流场对于科学合理地管理地下水系统的补给与排泄、恢复地下水位、修复相关生态环境问题具有重要支撑作用，已成为当前水资源保护和管理工作的重要支撑内容。
在地下水系统数值仿真中，往往将地下水含水层离散成具有指定行、列、层数的网格单元系统进行模拟，每个网格单元代表一定范围内的地下水含水层。目前，国际代表性地下水系统数值仿真模型如MODFLOW对大降深疏干条件下的地下水系统的模拟能力不强。大降深疏干条件下地下水系统模拟的主要特点在于会出现大量被疏干的地下水系统数值仿真网格单元(以下简称“网格单元”，图1)。在经典的MODFLOW-2005模型之中，疏干单元的出现会引发两方面问题，一是MODFLOW-2005将疏干单元处理为无效网格单元，地下水的流动会在疏干单元处中断，地下水的流动在疏干单元(含水层)处不连续，这与实际情况不相符合；同时对于疏干单元，也不再模拟作用于该疏干单元(含水层)之上的源汇项，因此当大量网格单元疏干时会导致地下水系统水量平衡、地下水水头模拟结果失真。尽管MODFLOW-2005可采用试算方法将符合预设湿润条件的疏干单元重新湿润使之再次成为有效单元，一定程度上可以减轻上述两方面问题带来的影响，但疏干单元的存在仍然在很大程度上限制了地下水系统数值仿真模型的应用面和准确性。为克服MODFLOW-2005对疏干单元的处理方式所引起的弊端，在MODFLOW的最新版本MODFLOW-NWT中采用了特殊的“上风格式”算法来计算网格单元间的渗流量，据此可保证地下水的流动不在疏干单元处中断，疏干单元仍然可以被处理为有效网格单元，作用于疏干单元之上的源汇项模拟不会被取消，但由于“上风格式”算法本身精度较低，因此也从另一个方面损害了地下水系统数值仿真模型的准确性。
发明内容
为解决上述技术问题，本发明提出了一种地下水大降深疏干下的全有效网格单元流场模拟方法，提升了地下水系统数值仿真在大降深疏干条件下的物理意义和准确性，拓宽了地下水系统数值仿真的应用面。
为实现上述目的，本发明提供了一种地下水大降深疏干下的全有效网格单元流场模拟方法，包括以下步骤：
对地下水系统进行网格化剖分，获取地下水系统数值仿真网格单元；
基于达西定律结合调和平均渗透系数和加权平均过水断面厚度，获取沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数；
基于达西定律结合调和平均越流系数，获取沿垂向上地下水系统数值仿真所述第一网格单元和第三网格单元之间的水力传导系数；
基于水量平衡原则，获取地下水系统数值仿真所述第一网格单元的贮水项；
基于所述沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数、所述沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数及所述地下水系统数值仿真所述第一网格单元的贮水项，结合所述第一网格单元的水量平衡关系，构建大降深疏干条件下描述地下水运动的模型，迭代获取地下水系统各数值仿真所述第一网格单元的水头，通过图像化处理后获取大降深疏干条件下的流场。
可选的，获取所述地下水系统数值仿真网格单元的方法包括：
将地下水含水层空间上离散为具有一定行数、列数、层数的数值仿真网格单元系统，每个地下水系统数值仿真网格单元代表一定范围内的地下水含水层。
可选的，所述调和平均渗透系数和所述加权平均过水断面厚度的获取方法包括：
根据网格单元水平向的渗透系数、行间距、列间距、水头和顶底板高程，通过计算获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的所述调和平均渗透系数和所述加权平均过水断面厚度。
可选的，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的所述调和平均渗透系数的方法包括：
其中，KRave为沿水平方向上第一网格单元和第二网格单元之间的所述调和平均渗透系数；KRi,j,k和KRi,j+1,k分别为第一网格单元和第二网格单元沿行方向上的渗透系数；ΔRj+1/2为第一网格单元和第二网格单元中心点之间的距离；ΔRj和ΔRj+1分别为第j列、第j+1列网格单元的行间距。
可选的，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的加权平均过水断面厚度的方法包括：
在含水层没有顶板高程限制条件时，所述第一网格单元和所述第二网格单元之间的平均过水断面厚度通过网格单元各自的过水断面厚度线性加权计算如下，
其中,HTi,j,k和HTi,j+1,k分别为第一网格单元和第二网格单元各自的过水断面厚度；hi,j,k和hi,j+1,k分别为这两个网格单元的水头；MaxB为这两个网格单元中相对较高的底板高程；HTave为第一网格单元和第二网格单元之间的加权平均过水断面厚度；
在含水层有顶板高程限制条件时，则网格单元间的加权平均过水断面厚度计算如下：
其中，Topi,j,k和Topi,j+1,k分别为第一网格单元和第二网格单元各自的顶板高程；Boti,j,k和Boti,j+1,k分别为第一网格单元和第二网格单元各自的底板高程。
可选的，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的水力传导系数的方法包括：
其中，CRi,j+1/2,k为沿行方向上第一网格单元和第二网格单元之间的水力传导系数；ΔCi为第i行网格单元的列间距；hup为第一网格单元和第二网格单元中水头相对较高的网格单元的水头。
可选的，获取沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数的方法及计算包括：
基于所述第一网格单元和所述第三网格单元之间的调和平均越流系数，及所述第一网格单元和所述第三网格单元的行间距、列间距，通过计算获取沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数；
CVi,j,k+1/2＝VCi,j,k+1/2·ΔRj·ΔCi
其中，CVi,j,k+1/2为沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数；VCi,j,k+1/2为第一网格单元和第三网格单元之间的调和平均越流系数。
可选的，获取地下水系统数值仿真所述第一网格单元的贮水项的方法包括：
SCIi,j,k＝ui,j,k·ΔRj·ΔCi
其中，SCIi,j,k为地下水系统数值仿真所述第一网格单元的贮水项，ui,j,k为贮水系数。
可选的，所述第一网格单元的水量平衡关系包括：
基于所述第一网格单元相邻的六个网格单元流入所述第一网格单元的水量与作用在所述第一网格单元之上的源汇项之和，去掉所述第一网格单元流出到所述六个网格单元的水量，获取所述第一网格单元的水量平衡关系。
可选的，构建大降深疏干条件下描述地下水运动的模型包括：
其中，下标(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为模拟时段内的水头变化量。
本发明技术效果：本发明公开了一种地下水大降深疏干下的全有效网格单元流场模拟方法，首先对地下水系统进行网格化剖分，得到地下水系统数值仿真网格单元；然后根据达西定律，结合调和平均渗透系数、加权平均过水断面厚度、调和平均越流系数计算地下水系统数值仿真网格单元间的水平向和垂向水力传导系数，以表征含水层中地下水侧向、垂向流动的能力；然后根据水量平衡原理计算网格单元的贮水项；最后在上述计算结果的基础上根据水量平衡原理构建大降深疏干条件下的描述地下水运动的差分方程组，求解得到大降深疏干条件下的模拟流场。和现有技术相比，地下水的流动在疏干单元处不会被人为中断，疏干单元仍被处理为有效网格单元，作用于疏干单元之上的源汇项模拟不会被取消，同时网格单元间渗流量算法拥有很高的计算精度，因此提升了地下水系统数值仿真在大降深疏干条件下的物理意义和准确性，拓宽了地下水系统数值仿真的应用面。
附图说明
构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：
图1为本发明实施例中网格单元被疏干的示意图；
图2为本发明实施例地下水大降深疏干下的全有效网格单元流场模拟方法的流程示意图；
图3为本发明实施例含水层的网格化空间离散技术示意图，其中(a)为含水层水平方向上(行方向、列方向)的网格化离散情况，(b)为含水层垂向上的网格化离散情况；
图4为本发明实施例沿水平方向上地下水系统数值仿真网格单元之间的水力传导系数计算流程图；
图5为本发明实施例沿垂向上地下水系统数值仿真网格单元之间的水力传导系数计算流程图；
图6为本发明实施例地下水系统数值仿真网格单元的贮水项计算流程图；
图7为本发明实施例任意一个网格单元与其周围六个相邻网格单元的示意图；
图8为本发明实施例使用Picard迭代法求解地下水运动差分方程组的计算流程图；
图9为本发明实施例模拟流场与解析流场的对比图。
具体实施方式
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。
如图1-9所示，本实施例中提供一种地下水大降深疏干下的全有效网格单元流场模拟方法，包括以下步骤：
S1、对地下水系统进行网格化剖分，得到地下水系统数值仿真网格单元；
S2、根据达西定律，结合调和平均渗透系数和加权平均过水断面厚度计算沿水平方向上地下水系统数值仿真网格单元之间的水力传导系数，以表征含水层中地下水侧向流动的能力。
S3、根据达西定律，结合调和平均越流系数计算沿垂向上地下水系统数值仿真网格单元之间的水力传导系数，以表征含水层中地下水垂向流动的能力。
S4、根据水量平衡原则，计算地下水系统数值仿真网格单元的贮水项。
S5、在S1、S2、S3、S4步骤基础上，结合源汇项条件，构建大降深疏干条件下描述地下水运动的差分方程组，采用Picard迭代法对该方程组进行迭代求解，得到地下水系统各数值仿真网格单元的水头，经图像化处理后即可得到大降深疏干条件下的流场。
在实施例中，模拟方法主要包括五部分内容，具体为，1，地下水系统的网格化剖分、2水平向水力传导系数计算、3垂向水力传导系数计算、4贮水项计算、5地下水运动差分方程组的构建与求解、模拟结果的图像化后处理，下面分别针对五部分内容进行详细的解释说明如图2。
一、地下水系统的网格化剖分
该部分内容对应于步骤S1，如图3所示，将地下水含水层空间上离散为具一定行数、列数、层数的数值仿真网格单元系统，每个地下水系统数值仿真网格单元代表一定范围内的地下水含水层。
二、水平向水力传导系数计算
该部分内容对应步骤S2，如图4所示，步骤S2具体为，①根据网格单元水平向(行方向、列方向)的渗透系数、行间距、列间距计算沿水平方向上(行方向、列方向)网格单元间的调和平均渗透系数；②根据网格单元的水头、顶底板高程、行间距、列间距计算沿水平方向(行方向、列方向)上网格单元间的加权平均过水断面厚度；③根据达西定律，结合单元间调和平均渗透系数、加权平均过水断面厚度计算结果计算沿水平方向上(行方向、列方向)网格单元间的水力传导系数，以表征含水层中地下水侧向流动的能力。沿行方向与沿列方向上地下水系统数值仿真网格单元间的调和平均渗透系数、平均过水断面厚度、水力传导系数的计算方法完全一致，以行方向为例上述变量分别计算如下：
1、对任意一个网格单元(i,j,k)和其行方向上直接相邻的网格单元(i,j+1,k)，利用下式计算沿行方向上网格单元间的调和平均渗透系数，
其中,KRi,j,k和KRi,j+1,k分别为网格单元(i,j,k)和(i,j+1,k)沿行方向上的渗透系数(LT-1)；ΔRj+1/2为这两个网格单元中心点之间的距离(L)；ΔRj和ΔRj+1分别为第j列、第j+1列网格单元的行间距(L)。
2、在含水层没有顶板高程限制条件时，网格单元(i,j,k)和(i,j+1,k)之间的平均过水断面厚度HTave(L)通过网格单元各自的过水断面厚度线性加权计算如下，
其中，HTi,j,k和HTi,j+1,k分别为网格单元(i,j,k)和(i,j+1,k)各自的过水断面厚度(L)；hi,j,k和hi,j+1,k分别为这两个网格单元的水头(L)；MaxB为这两个网格单元中相对较高的底板高程(L)。
在含水层有顶板高程限制条件时，则网格单元间的加权平均过水断面厚度计算如下：
其中，Topi,j,k和Topi,j+1,k分别为网格单元(i,j,k)和(i,j+1,k)各自的顶板高程(L)；Boti,j,k和Boti,j+1,k分别为这两个网格单元各自的底板高程(L)；
3、根据达西定律，在调和平均渗透系数、加权平均过水断面厚度计算结果的基础上，利用下式计算沿行方向上网格单元(i,j,k)和(i,j+1,k)之间的水力传导系数，以表征含水层中地下水侧向流动的能力，
其中，CRi,j+1/2,k为沿行方向上网格单元(i,j,k)和(i,j+1,k)之间的水力传导系数(L2T-1)；ΔCi为第i行网格单元的列间距(L)；hup为这两个网格单元中水头相对较高的网格单元的水头(L)。
三、垂向水力传导系数计算
该部分内容对应步骤S3,如图5所示，步骤S3具体为，①判断是否已知上下层网格单元之间的调和平均越流系数；②若已知上下层网格单元之间的调和平均越流系数则根据达西定律结合网格单元的行间距、列间距直接计算上下层网格单元之间的垂向水力传导系数；若上下层网格单元之间的调和平均越流系数是未知的，则判断上下层网格单元之间是否存在低渗透性介质；③若不存在低渗透性介质，则根据网格单元的垂向渗透系数、水头、顶底板高程先计算上下层网格单元间的调和平均越流系数，再结合网格单元的行间距、列间距计算上下层网格单元之间的垂向水力传导系数；若存在低渗透性介质，则根据网格单元的垂向渗透系数、水头、顶底板高程、上下层网格单元间低渗透性介质的厚度和渗透系数先计算上下层网格单元间的调和平均越流系数，再结合网格单元的行间距、列间距计算上下层网格单元之间的垂向水力传导系数。上下层网格单元间的调和平均越流系数、垂向水力传导系数分别计算如下：
1、若已知上下层网格单元之间的调和平均越流系数，根据下式计算上下层网格单元之间的垂向水力传导系数
CVi,j,k+1/2＝VCi,j,k+1/2·ΔRj·ΔCi
其中，CVi,j,k+1/2为任意一个上层单元(i,j,k)和其直接相邻的下层网格单元(i,j,k+1)之间的垂向水力传导系数(L2T-1)；VCi,j,k+1/2为上层网格单元(i,j,k)和下层网格单元(i,j,k+1)之间的调和平均越流系数(T-1)；
2、若上下层网格单元之间的调和平均越流系数是未知的，则先计算沿垂向上上下层网格单元间的调和平均越流系数，
HSati,j,k＝min[max(0,hi,j,k-Boti,j,k),Topi,j,k-Boti,j,k]
TKi,j,k+1＝Topi,j,k+1-Boti,j,k+1
其中，HSati,j,k为上层网格单元(i,j,k)的含水层饱和厚度(L)，当上层网格单元承压时，该值等于该网格单元的含水层厚度(L)；TKi,j,k+1为下层网格单元(i,j,k+1)的含水层厚度(L)；TKCB为上、下层网格单元之间低渗透性介质的厚度(L)；VKi,j,k、VKi,j,k+1及VKB分别为上层网格单元、下层网格单元及低渗透性介质的垂向渗透系数(LT-1)；TOPi,j,k+1和Boti,j,k+1分别为下层网格单元(i,j,k+1)的顶、底板高程值(L)。
若上下层网格单元间不存在低渗透性介质，则忽略上式中的TKCB和VKB项。
根据上、下层网格单元之间的调和平均越流系数可进一步计算上、下层网格单元之间的垂向水力传导系数，以表征含水层中地下水垂向流动的能力。
四、贮水项计算
该部分内容对应步骤S4,如图6所示，步骤S4具体为，对任意一个网格单元(i,j,k)采用不同的贮水系数分别计算其处于疏干、湿润-非承压、湿润-承压三种状态下的贮水项，以表征网格单元(含水层)水头上升或下降一个单位时，网格单元(含水层)处于不同状态时存入的或释放的水量，
SCIi,j,k＝ui,j,k·ΔRj·ΔCi
其中，SCIi,j,k为网格单元(i,j,k)的贮水项(L2)，ui,j,k为贮水系数。
若网格单元(i,j,k)的模拟水头不高于其底板高程，即该单元处于疏干状态时，由水量平衡原理可知，该第一网格单元虽有计算意义上的水头，但该单元贮水量为0，此时ui,j,k＝0；
若网格单元(i,j,k)的模拟水头高于其底板高程低于其顶板高程，即该单元处于湿润-非承压状态时，贮水系数即为给水度，即ui,j,k＝SY，SY表示给水度。
若网格单元(i,j,k)的模拟水头不低于其顶板高程，即该单元处于承压状态时，贮水系数大于0，即ui,j,k＞0。
五、地下水运动差分方程组的构建与求解、模拟结果的图像化后处理
该部分内容对应步骤S5，步骤S5具体为，对任一网格单元(i,j,k)，在一定时间内从其相邻的6个网格单元(上、下、前、后、左、右)流入第一网格单元(i,j,k)的水量与作用于第一网格单元(i,j,k)之上的源汇项之和，减去从第一网格单元(i,j,k)流出到六个相邻网格单元的水量应该等于网格单元(i,j,k)的贮水量的变化(图7)，据此水量平衡关系结合步骤S2、S3、S4的计算结果建立有关网格单元(i,j,k)处地下水运动的差分方程，
其中，下标(i,j,k-1)，(i,j,k+1)，(i+1,j,k)，(i-1,j,k)，(i,j-1,k)，(i,j+1,k)，分别代表网格单元(i,j,k)的上、下、前、后、左、右侧的相邻网格单元(图7)；CR、CC、CV分别为沿行、列、层方向上网格单元之间的水力传导系数(L2T-1)；h为水头(L)；上标n代表本次迭代，上标n-1代表上次迭代；Pi,j,k代表作用于网格单元(i,j,k)之上的与水头有关的源汇项相关系数(L2T-1)；Qi,j,k代表作用于网格单元(i,j,k)之上的流量源汇项(L3T-1)；SCIi,j,k为网格单元(i,j,k)的贮水项(L2)；Δt代表模拟时段的时间长度，Δh代表模拟时段内的水头变化量(L)。
对地下水系统所涉及的所有网格单元逐个列出以上差分方程，形成刻画整个地下水系统中地下水运动的差分方程组，求解即得模拟时段末整个地下水系统的模拟水头分布。
对地下水运动差分方程组采用Picard迭代法求解，已知地下水系统各数值仿真网格单元的初始水头值h0，需模拟t时间后的水头值ht,Picard迭代法求解的步骤可描述如图8，①根据h0分别计算计算沿水平向和垂向上网格单元之间的水力传导系数CR、CC、CV；②利用0值贮水系数、给水度、大于0的贮水系数分别计算网格单元处于疏干、湿润-非承压、湿润-承压三种状态下的贮水项；③结合源汇项构建描述地下水运动的差分方程组，并求解，得到第一次迭代计算求解的水头h1；④判断所有网格单元上第一次迭代求解的水头值和初始水头值相差(|h1-h0|)是否小于事先规定的阈值hclose；⑤如果是，表示模拟收敛，此时模拟结束，ht＝h1，如果不是，用h1取代h0更新地下水运动差分方程的相关系数项步骤①③，再次求解得到第二次迭代求解的水头h2，若在所有网格上|h2-h1|仍不小于hclose，则继续重复步骤用h2取代h1，计算h3,h4…hn,直到满足在所有网格单元上本次迭代和上次迭代求解的水头值|hn-hn-1|＜hclose，此时认为模拟收敛，不再进行迭代求解，模拟结束，即ht＝hn。
最后，采用Origin、Surfer等类似绘图软件对求解的各网格单元水头值进行图像化处理即可获得大降深疏干条件下的模拟地下水流场。
本实施例情况：
某地含水层总面积25km2，含水层厚度为50米，含水层东侧水头为50m，西侧水头10m,含水层东西侧存在40m的大降深，地下水自东向西流出含水层，要求模拟稳定状态下该含水层的水头值；将该含水层空间上离散为20行、20列、50层的网格单元系统(图2)。作为对比，也使用MODFLOW-2005和MODFLOW-NWT模拟本实施例的含水层。
达到的应用效果：
受大降深影响，本实施例有超过30％的网格在模拟结束后仍处于疏干状态。图9为本发明的模拟水头误差与MODFLOW-2005、MODFLOW-NWT的模拟水头误差对比，可以看出本发明的模拟水头与解析水头几乎完全一致，MODFLOW-2005和MODFLOW-NWT的模拟水头误差明显高于本发明，最大分别超过1m和0.6m。流量模拟结果如表1所示，本发明的模拟流量与解析流量几乎完全一致，MODFLOW-2005和MODFLOW-NWT的模拟流量与解析流量有明显偏差。
表1
以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。
Claims (10)
1.一种地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，包括以下步骤：
对地下水系统进行网格化剖分，获取地下水系统数值仿真网格单元；
基于达西定律结合调和平均渗透系数和加权平均过水断面厚度，获取沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数；
基于达西定律结合调和平均越流系数，获取沿垂向上地下水系统数值仿真所述第一网格单元和第三网格单元之间的水力传导系数；
基于水量平衡原则，获取地下水系统数值仿真所述第一网格单元的贮水项；
基于所述沿水平方向上地下水系统数值仿真第一网格单元和第二网格单元之间的水力传导系数、所述沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数及所述地下水系统数值仿真所述第一网格单元的贮水项，结合所述第一网格单元的水量平衡关系，构建大降深疏干条件下描述地下水运动的模型，迭代获取地下水系统各数值仿真所述第一网格单元的水头，通过图像化处理后获取大降深疏干条件下的流场。
2.如权利要求1所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，获取所述地下水系统数值仿真网格单元的方法包括：
将地下水含水层空间上离散为具有一定行数、列数、层数的数值仿真网格单元系统，每个地下水系统数值仿真网格单元代表一定范围内的地下水含水层。
3.如权利要求1所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，所述调和平均渗透系数和所述加权平均过水断面厚度的获取方法包括：
根据网格单元水平向的渗透系数、行间距、列间距、水头和顶底板高程，通过计算获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的所述调和平均渗透系数和所述加权平均过水断面厚度。
4.如权利要求3所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的所述调和平均渗透系数的方法包括：
其中，KRave为沿水平方向上第一网格单元和第二网格单元之间的所述调和平均渗透系数；KRi,j,k和KRi,j+1,k分别为第一网格单元和第二网格单元沿行方向上的渗透系数；ΔRj+1/2为第一网格单元和第二网格单元中心点之间的距离；ΔRj和ΔRj+1分别为第j列、第j+1列网格单元的行间距。
5.如权利要求3所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的加权平均过水断面厚度的方法包括：
在含水层没有顶板高程限制条件时，所述第一网格单元和所述第二网格单元之间的平均过水断面厚度通过网格单元各自的过水断面厚度线性加权计算如下，
其中,HTi,j,k和HTi,j+1,k分别为第一网格单元和第二网格单元各自的过水断面厚度；hi,j,k和hi,j+1,k分别为这两个网格单元的水头；MaxB为这两个网格单元中相对较高的底板高程；HTave为第一网格单元和第二网格单元之间的加权平均过水断面厚度；
在含水层有顶板高程限制条件时，则网格单元间的加权平均过水断面厚度计算如下：
其中，Topi,j,k和Topi,j+1,k分别为第一网格单元和第二网格单元各自的顶板高程；Boti,j,k和Boti,j+1,k分别为第一网格单元和第二网格单元各自的底板高程。
6.如权利要求3所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，获取沿水平方向上地下水系统数值仿真所述第一网格单元和所述第二网格单元之间的水力传导系数的方法包括：
其中，CRi,j+1/2,k为沿行方向上第一网格单元和第二网格单元之间的水力传导系数；ΔCi为第i行网格单元的列间距；hup为第一网格单元和第二网格单元中水头相对较高的网格单元的水头。
7.如权利要求1所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，获取沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数的方法及计算包括：
基于所述第一网格单元和所述第三网格单元之间的调和平均越流系数，及所述第一网格单元和所述第三网格单元的行间距、列间距，通过计算获取沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数；
CVi,j,k+1/2＝VCi,j,k+1/2·ΔRj·ΔCi
其中，CVi,j,k+1/2为沿垂向上地下水系统数值仿真所述第一网格单元和所述第三网格单元之间的水力传导系数；VCi,j,k+1/2为第一网格单元和第三网格单元之间的调和平均越流系数。
8.如权利要求1所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，获取地下水系统数值仿真所述第一网格单元的贮水项的方法包括：
SCIi,j,k＝ui,j,k·ΔRj·ΔCi
其中，SCIi,j,k为地下水系统数值仿真所述第一网格单元的贮水项，ui,j,k为贮水系数。
9.如权利要求1所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，所述第一网格单元的水量平衡关系包括：
基于所述第一网格单元相邻的六个网格单元流入所述第一网格单元的水量与作用在所述第一网格单元之上的源汇项之和，去掉所述第一网格单元流出到所述六个网格单元的水量，获取所述第一网格单元的水量平衡关系。
10.如权利要求1所述的地下水大降深疏干下的全有效网格单元流场模拟方法，其特征在于，构建大降深疏干条件下描述地下水运动的模型包括：
其中，下标(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为模拟时段内的水头变化量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310582904.9A CN116992783B (zh) | 2023-05-23 | 2023-05-23 | 一种地下水大降深疏干下的全有效网格单元流场模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310582904.9A CN116992783B (zh) | 2023-05-23 | 2023-05-23 | 一种地下水大降深疏干下的全有效网格单元流场模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116992783A true CN116992783A (zh) | 2023-11-03 |
| CN116992783B CN116992783B (zh) | 2024-06-14 |
Family
ID=88525509
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310582904.9A Active CN116992783B (zh) | 2023-05-23 | 2023-05-23 | 一种地下水大降深疏干下的全有效网格单元流场模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116992783B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102156779A (zh) * | 2011-04-13 | 2011-08-17 | 北京石油化工学院 | 地下水流仿真与预测分析方法 |
| US20130333881A1 (en) * | 2012-06-14 | 2013-12-19 | Besst, Inc. | Selective extraction of fluids from subsurface wells |
| CN105022913A (zh) * | 2015-06-01 | 2015-11-04 | 中国水利水电科学研究院 | 一种降雨入渗补给地下水临界埋深计算方法 |
| JP2016017773A (ja) * | 2014-07-04 | 2016-02-01 | 国立大学法人信州大学 | 熱応答試験および揚水試験の解析方法および解析プログラム |
| CN108763797A (zh) * | 2018-06-04 | 2018-11-06 | 中国水利水电科学研究院 | 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 |
| CN113536644A (zh) * | 2021-07-30 | 2021-10-22 | 江苏南京地质工程勘察院 | 悬挂式止水帷幕深基坑降水方案模拟优化方法 |
| US20230115283A1 (en) * | 2021-10-08 | 2023-04-13 | Hohai University | Three-level grid multi-scale finite element method for simulating groundwater flow in heterogeneous aquifers |
-
2023
- 2023-05-23 CN CN202310582904.9A patent/CN116992783B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102156779A (zh) * | 2011-04-13 | 2011-08-17 | 北京石油化工学院 | 地下水流仿真与预测分析方法 |
| US20130333881A1 (en) * | 2012-06-14 | 2013-12-19 | Besst, Inc. | Selective extraction of fluids from subsurface wells |
| JP2016017773A (ja) * | 2014-07-04 | 2016-02-01 | 国立大学法人信州大学 | 熱応答試験および揚水試験の解析方法および解析プログラム |
| CN105022913A (zh) * | 2015-06-01 | 2015-11-04 | 中国水利水电科学研究院 | 一种降雨入渗补给地下水临界埋深计算方法 |
| CN108763797A (zh) * | 2018-06-04 | 2018-11-06 | 中国水利水电科学研究院 | 一种基于地下水模型的湖泊与地下水稳定流作用模拟方法 |
| CN113536644A (zh) * | 2021-07-30 | 2021-10-22 | 江苏南京地质工程勘察院 | 悬挂式止水帷幕深基坑降水方案模拟优化方法 |
| US20230115283A1 (en) * | 2021-10-08 | 2023-04-13 | Hohai University | Three-level grid multi-scale finite element method for simulating groundwater flow in heterogeneous aquifers |
Non-Patent Citations (3)
| Title |
|---|
| 吴昌友;李中才;王福林;: "三江平原挠力河流域地下水流数学模型的建立及仿真", 水动力学研究与进展A辑, no. 03 * |
| 孟宪萌;邵骏煜;尹茂生;刘登峰;: "越流系统水流运动规律研究综述", 水利水电科技进展, no. 04 * |
| 赵文凤;谢一凡;吴吉春;: "一种模拟节点达西渗透流速的双重网格多尺度有限单元法", 岩土工程学报, no. 08, 15 August 2020 (2020-08-15) * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116992783B (zh) | 2024-06-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109492299B (zh) | 基于swmm与modflow耦合的水资源模拟方法 | |
| CN107832931B (zh) | 一种平原水网地区内涝风险的模块化分析方法 | |
| Sedki et al. | Simulation-optimization modeling for sustainable groundwater development: a Moroccan coastal aquifer case study | |
| Shrestha et al. | Mapping groundwater resiliency under climate change scenarios: a case study of Kathmandu Valley, Nepal | |
| Haregeweyn et al. | Assessing the performance of a spatially distributed soil erosion and sediment delivery model (WATEM/SEDEM) in Northern Ethiopia | |
| Ahmed et al. | Groundwater flow modelling of Yamuna-Krishni interstream, a part of central Ganga Plain Uttar Pradesh | |
| Bhandari et al. | 2D unsteady flow routing and flood inundation mapping for lower region of Brazos River watershed | |
| CN106971034B (zh) | 一种基于集合卡尔曼滤波的对无资料地区径流量推求方法 | |
| CN111062125B (zh) | 海绵型综合管廊水文效应评估方法 | |
| Alvarez et al. | Groundwater flow model, recharge estimation and sustainability in an arid region of Patagonia, Argentina | |
| Abbas et al. | Improving river flow simulation using a coupled surface-groundwater model for integrated water resources management | |
| CN112149871A (zh) | 一种基于gis空间统计与随机模拟相结合的污染物点源解析方法 | |
| CN112907047A (zh) | 一种海绵城市绩效评估系统 | |
| Efstratiadis et al. | HYDROGEIOS: a semi-distributed GIS-based hydrological model for modified river basins | |
| CN114004141A (zh) | 一种河流流域地下水超采区预测决策平台 | |
| Chang et al. | Assessing the impact of climate variability and human activity to streamflow variation | |
| Bemporad et al. | A distributed approach for sediment yield evaluation in Alpine regions | |
| CN116992783B (zh) | 一种地下水大降深疏干下的全有效网格单元流场模拟方法 | |
| He et al. | Regional groundwater prediction model using automatic parameter calibration SCE method for a coastal plain of Seto Inland Sea | |
| Tîrnovan et al. | Predicting the potential index of major floods production in the Suha river basin (Suha Bucovineana) | |
| Wijaya et al. | Flood Mapping Using HEC-RAS and HEC-HMS: A Case Study of Upper Citarum River at Dayeuhkolot District, Bandung Regency, West Java | |
| Charbonneau et al. | Wetland Modeling in PCSWMM: Exploring Options to Define Wetland Features and Incorporate Groundwater Exchanges | |
| Ghorbani et al. | Interceptor Drainage Modelling to Manage High Groundwater Table on the Abyek Plain, Iran | |
| Manaouch et al. | Integrating WaTEM/SEDEM model and GIS-based FAHP Method for Identifying Ecological Rainwater Harvesting Sites in Ziz upper watershed, SE Morocco | |
| CN110717247B (zh) | 一种考虑梯田影响的坡面产汇流过程模拟方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information | ||
| CB03 | Change of inventor or designer information |
Inventor after: Lu Chuiyu Inventor after: Wu Weichen Inventor after: Li Xinwang Inventor after: Qin Tao Inventor after: Wu Chu Inventor after: Yan Lingjia Inventor after: Lu Wen Inventor after: He Xin Inventor after: Sun Qingyan Inventor after: Liu Miao Inventor after: Wu Zhenjiang Inventor after: Han Shangqi Inventor before: Lu Chuiyu Inventor before: Wu Weichen Inventor before: Liu Xinwang Inventor before: Qin Tao Inventor before: Wu Chu Inventor before: Yan Lingjia Inventor before: Lu Wen Inventor before: He Xin Inventor before: Sun Qingyan Inventor before: Liu Miao Inventor before: Wu Zhenjiang Inventor before: Han Shangqi |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |