CN116522818A - 一种大坡度地形条件下的干旱区潜水位模拟方法 - Google Patents

一种大坡度地形条件下的干旱区潜水位模拟方法 Download PDF

Info

Publication number
CN116522818A
CN116522818A CN202310512745.5A CN202310512745A CN116522818A CN 116522818 A CN116522818 A CN 116522818A CN 202310512745 A CN202310512745 A CN 202310512745A CN 116522818 A CN116522818 A CN 116522818A
Authority
CN
China
Prior art keywords
grid
iteration
diving
simulation
adjacent
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
Application number
CN202310512745.5A
Other languages
English (en)
Other versions
CN116522818B (zh
Inventor
陆垂裕
严聆嘉
吴初
陆文
何鑫
孙青言
刘淼
吴镇江
韩尚麒
吴委尘
张召召
秦韬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Water Resources and Hydropower Research
Original Assignee
China Institute of Water Resources and Hydropower Research
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Institute of Water Resources and Hydropower Research filed Critical China Institute of Water Resources and Hydropower Research
Priority to CN202310512745.5A priority Critical patent/CN116522818B/zh
Publication of CN116522818A publication Critical patent/CN116522818A/zh
Application granted granted Critical
Publication of CN116522818B publication Critical patent/CN116522818B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Fluid Mechanics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出了一种大坡度地形条件下的干旱区潜水位模拟方法，包括：S1、对潜水含水层进行网格化空间离散，获取潜水含水层模拟网格单元；S2、构建潜水数值仿真矩阵方程，求解得到本次迭代下各所述模拟网格单元的模拟潜水位；S3、判断本次迭代下各所述模拟网格单元的模拟潜水位是否满足预设精度阈值，若满足则结束模拟，否则返回S2进行下一次迭代。本发明提升了大坡度地形条件下干旱区潜水位数值仿真的成功率和准确性，拓宽了地下水数值仿真的应用面。

Description

一种大坡度地形条件下的干旱区潜水位模拟方法
技术领域
本发明属于地下水数值仿真技术领域，尤其涉及一种大坡度地形条件下的干旱区潜水位模拟方法。
背景技术
潜水位是指潜水含水层的水位，也即埋藏在地面以下第一个稳定隔水层之上具有自由重力水面的含水层的水位。干旱区潜水位直接影响干旱区植被的生长发育、盖度、年龄结构、种群构成与演替和生物多样性，准确模拟干旱区潜水位对于维持干旱区生态系统的稳定具有重要意义。我国西北部干旱区常伴随着大坡度地形，然而目前已有地下水数值仿真技术对大坡度地形条件下的干旱区潜水位模拟存在困难。
目前，国际最具代表性的地下水数值仿真模型是美国地调局开发的MODFLOW-2005模型。MODFLOW-2005将潜水含水层离散成一系列网格单元，每个网格单元代表一定范围内的潜水含水层；根据达西定律可计算相邻网格单元(含水层各部分)之间的地下水渗流量，据此可构建每个网格单元(含水层各部分)的水量平衡方程，通过联立各网格单元的水量平衡方程可形成描述整个潜水含水层中地下水运动的方程组(即潜水数值仿真矩阵方程)，求解可得各网格单元在某一时刻的潜水位，从而达到反演或预测潜水含水层系统中各处潜水位的目的。然而，以上模拟方法应用于大坡度地形条件下的干旱区潜水位模拟时会出现两方面问题，一是在大坡度地形条件下，水平方向上相邻网格单元之间很容易出现具有不完全水力联系或完全不具有水力联系的情况，此时传统的达西定律实际上并不适用，模型物理机制存在不足。二是受干旱区降水入渗补给匮乏、水文地质条件、人类活动等因素综合影响，干旱区潜水含水层各部分的水层厚度(潜水位减含水层底板高程的差值)急剧变化，使得潜水数值仿真矩阵方程的非线性程度较高，此时不合理地应用传统达西定律极易导致潜水数值仿真矩阵方程难以求解，从而导致模型模拟失败。
发明内容
为解决上述技术问题，本发明提供一种大坡度地形条件下的干旱区潜水位模拟方法，提升了大坡度地形条件下干旱区潜水位数值仿真的成功率和准确性，拓宽了地下水数值仿真的应用面。
为实现上述目的，本发明提供了一种大坡度地形条件下的干旱区潜水位模拟方法，包括：
S1、对潜水含水层进行网格化空间离散，获取潜水含水层模拟网格单元；
S2、构建潜水数值仿真矩阵方程，求解本次迭代下各所述模拟网格单元的模拟潜水位；
S3、判断本次迭代下各所述模拟网格单元的模拟潜水位是否满足预设精度阈值，若满足则结束模拟，否则返回S2进行下一次迭代。
可选地，对潜水含水层进行网格化空间离散包括：对潜水含水层从行方向、列方向和层方向进行网格化空间离散。
可选地，构建所述潜水数值仿真矩阵方程包括：
基于水量平衡关系，构建描述所述模拟网格单元处潜水运动的水量平衡方程；其中，所述水量平衡关系为：在模拟时段内自周围六个相邻网格单元流入模拟网格单元(i,j,k)的潜水量与作用于模拟网格单元(i,j,k)之上的源汇项之和等于模拟网格单元(i,j,k)贮水量的变化；
联立潜水含水层中全部所述模拟网格单元的水量平衡方程，获取所述潜水数值仿真矩阵方程。
可选地，所述水量平衡方程为：
其中，qi,j-1/2,k、qi,j+1/2,k、qi-1/2,j,k、qi+1/2,j,k、qi,j,k-1/2、qi,j,k+1/2分别表示网格单元(i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k)、(i,j,k-1)、(i,j,k+1)和(i,j,k)之间的渗流量，渗流量为正值表示潜水流向网格单元(i,j,k)，渗流量为负值表示潜水流出网格单元(i,j,k)；QSi,j,k为作用于网格单元(i,j,k)上的源汇项；表示网格单元(i,j,k)贮水量的变化率，(i,j,k)表示处于第i行、第j列、第k层的模拟网格单元。
可选地，所述潜水数值仿真矩阵方程为：
[A]{h}＝{q}
其中，[A]为系数矩阵；{h}为待求解的潜水位向量；{q}为右端项向量，代表水量平衡方程中的常数项和已知项。
可选地，求解本次迭代下各所述模拟网格单元的模拟潜水位包括：
在水平方向上，基于预设方程组，对模拟网格单元(i,j,k)及任一相邻网格单元((i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k))之间的水平向渗流量进行计算；其中，(i,j-1,k)、(i,j+1,k)分别为沿行方向上左、右侧相邻的网格单元；(i+1,j,k)、(i-1,j,k)分别为沿列方向上前、后侧相邻的网格单元；
在层方向上，基于预设方程组，对模拟网格单元(i,j,k)及任一相邻的网格单元((i,j,k-1)、(i,j,k+1))之间的垂向渗流量进行计算；其中，其中(i,j,k-1)、(i,j,k+1)分别为沿层方向上上、下侧相邻的网格单元；
将水平方向上各所述模拟网格单元和对应相邻网格单元之间的渗流量、层方向上各所述模拟网格单元和对应相邻网格单元的之间的渗流量、各所述模拟网格单元的源汇项、各所述模拟网格单元贮水量的变化率代入潜水数值仿真矩阵方程求解获取本次迭代下各所述模拟网格单元的模拟潜水位。
可选地，在水平方向上，基于预设方程组，对模拟网格单元(i,j,k)及任一相邻的网格单元((i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k))之间的渗流量进行计算包括：
用分别表示网格单元(i,j,k)和该相邻网格单元上次迭代求解的潜水位，用/>分别表示本次迭代待求解的潜水位，用Boti,j,k、Botadjoin分别表示底板高程，用C表示网格单元(i,j,k)和该相邻网格单元之间的水力传导系数；
(1)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系，本次迭代中渗流量符合传统达西定律，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(2)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自网格单元(i,j,k)流向该相邻网格单元，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(3)当Boti,j,k≥Botadjoin时，若则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自该相邻网格单元流向网格单元(i,j,k)，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(4)当Boti,j,k≥Botadjoin时，若且/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系，单元间无渗流量：
q＝0
(5)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系，本次迭代中渗流量符合传统达西定律，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(6)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自该相邻网格单元流向网格单元(i,j,k)，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(7)当Botadjoin＞Boti,j,k时，若则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自网格单元(i,j,k)流向该相邻网格单元，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(8)当Botadjoin＞Boti,j,k时，若且/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系，单元间无渗流量：
q＝0。
可选地，本次迭代中网格单元(i,j,k)和沿层方向上的相邻网格单元(i,j,k-1)、(i,j,k+1)之间的垂向渗流量分别为：
其中，Ci,j,k-1/2、Ci,j,k+1/2分别表示网格单元(i,j,k)和(i,j,k-1)、(i,j,k+1)之间的水力传导系数，分别表示网格单元(i,j,k-1)和(i,j,k+1)本次迭代待求解的潜水位。
可选地，本次迭代作用于所述模拟网格单元之上的源汇项为：
其中，Pi,j,k代表作用于网格单元(i,j,k)之上的与本次迭代待求解潜水位有关的源汇项相关系数；Qi,j,k代表作用于网格单元(i,j,k)之上的流量源汇项；
本次迭代计算的所述模拟网格单元贮水量的变化率为：
其中，SSi,j,k代表网格单元(i,j,k)的贮水率，Vi,j,k代表网格单元(i,j,k)的体积，代表本次迭代计算的网格单元(i,j,k)的潜水位变化率。
可选地，判断本次迭代模拟的潜水位是否满足预设精度阈值包括：
将本次迭代模拟的潜水位与上次迭代模拟的潜水位进行比较，当所有所述模拟网格单元上本次迭代模拟的潜水位与上次迭代模拟的潜水位相差都小于预先规定的精度阈值时，认为模拟收敛，即潜水数值仿真矩阵方程求解完成，不再进行迭代。
与现有技术相比，本发明具有如下优点和技术效果：
本发明在求解潜水数值仿真矩阵方程时，利用修正后的达西定律确定水平方向上具有不完全水力联系或完全不具有水力联系的相邻网格单元之间的渗流量，提高了大坡度地形条件下干旱区潜水位数值仿真的成功率和准确性。
附图说明
构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：
图1为本发明实施例的一种大坡度地形条件下的干旱区潜水位模拟方法流程示意图；
图2为本发明实施例的基于矩形网格的潜水含水层空间离散技术示意图；
图3为本发明实施例的任意一个网格单元与其周围六个相邻网格单元的示意图；
图4为本发明实施例的网格单元(i,j,k)的底板高程不低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系的情况示意图；
图5为本发明实施例的网格单元(i,j,k)的底板高程不低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间第一种具有不完全水力联系的情况示意图；
图6为本发明实施例的网格单元(i,j,k)的底板高程不低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间第二种具有不完全水力联系的情况示意图；
图7为本发明实施例的网格单元(i,j,k)的底板高程不低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间完全不具有的水力联系的情况示意图；
图8为本发明实施例的网格单元(i,j,k)的底板高程低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系的情况示意图；
图9为本发明实施例的网格单元(i,j,k)的底板高程低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间第一种具有不完全水力联系的情况示意图；
图10为本发明实施例的网格单元(i,j,k)的底板高程低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间第二种具有不完全水力联系的情况示意图；
图11为本发明实施例的网格单元(i,j,k)的底板高程低于水平向相邻网格单元的底板高程时，网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系的情况示意图；
图12为本发明实施例的大坡度地形立体示意图；
图13为本发明实施例的模拟潜水位图。
具体实施方式
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。
本发明实施例提供了一种大坡度地形条件下的干旱区潜水位模拟方法，如图1所示，包括以下步骤S1-S3。
S1、对潜水含水层进行网格化空间离散，获取潜水含水层模拟网格单元；
进一步地，步骤S1具体为，采用矩形网格将潜水含水层空间上离散为具有一定行数、列数、层数的模拟网格单元系统(图2)；
S2、构建潜水数值仿真矩阵方程，求解得到本次迭代下各模拟网格单元的模拟潜水位；
进一步地，步骤S2具体为，不失一般的，对任意一个处于第i行、第j列、第k层的模拟网格单元(i,j,k)，其周围至多存在六个相邻网格单元(图3)，沿行方向上左、右侧相邻的网格单元分别为(i,j-1,k)、(i,j+1,k)，沿列方向上前、后侧相邻的网格单元分别为(i+1,j,k)、(i-1,j,k)，沿层方向上上、下侧相邻的网格单元分别为(i,j,k-1)、(i,j,k+1)，按照“在模拟时段内自周围六个相邻网格单元流入网格单元(i,j,k)的潜水量(可正可负)与作用于网格单元(i,j,k)之上的源汇项之和等于网格单元(i,j,k)贮水量的变化”这一水量平衡关系，建立如下描述网格单元(i,j,k)处潜水运动的水量平衡方程，
其中，qi,j-1/2,k(L3T-1)、qi,j+1/2,k(L3T-1)、qi-1/2,j,k(L3T-1)、qi+1/2,j,k(L3T-1)、qi,j,k-1/2(L3T-1)、qi,j,k+1/2(L3T-1)分别代表网格单元(i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k)、(i,j,k-1)、(i,j,k+1)和(i,j,k)之间的渗流量，正值代表潜水流向网格单元(i,j,k)，负值代表潜水流出网格单元(i,j,k)；QSi,j,k为作用于网格单元(i,j,k)上的源汇项(L3T-1)；代表网格单元(i,j,k)贮水量的变化率(L3T-1)。
进一步地，式(1)中各项都与待求解的潜水位有关，联立潜水含水层系统所有模拟网格单元的式(1)，整理成矩阵方程形式即得潜水数值仿真矩阵方程，
[A]{h}＝{q} (2)
其中，[A]为系数矩阵(L2T-1)；{h}为待求解的潜水位向量(L)，{q}为右端项向量(L3T-1)，代表式(1)中的常数项和已知项。
进一步地，潜水数值仿真矩阵方程的系数矩阵[A]与待求解的潜水位有关,因此矩阵方程无法直接求解，而是需要反复迭代计算直到求解的潜水位满足事先规定的精度要求。迭代求解的步骤为，若已知模拟时段初各模拟网格单元的初始潜水位组成的向量为h0，需要求解模拟时段末各模拟网格单元的潜水位，则首先利用h0构建系数矩阵[A]和右端项向量{q}，求解潜水数值仿真矩阵方程可得第一次迭代求解的各模拟网格单元的潜水位h1，之后再利用h1更新系数矩阵[A]和右端项向量{q}，再次求解潜水数值仿真矩阵方程可得第二次迭代求解的各模拟网格单元的潜水位h2，如此循环往复不断迭代，直到求解的潜水位满足事先规定的精度要求为止，最后一次迭代求解的潜水位就是模拟时段末的潜水位。
进一步地，以网格单元(i,j,k)的水量平衡方程(式1)为例说明潜水数值仿真矩阵方程各项计算方式。在水平方向上(行方向、列方向)，对网格单元(i,j,k)及其任一相邻的网格单元((i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k))，若这两个网格单元上次迭代求解的潜水位都高于这两个网格单元之中相对较高的底板高程，则认为本次迭代中这两个相邻网格单元之间具有完全水力联系；若这两个网格单元上次迭代求解的潜水位只有一个高于这两个网格单元之中相对较高的底板高程，则认为本次迭代中这两个相邻网格单元之间具有不完全的水力联系；若这两个网格单元上次迭代求解的潜水位都不高于这两个网格单元之中相对较高的底板高程，则认为本次迭代中这两个相邻网格单元之间完全不具有水力联系。不失一般性，考虑网格单元(i,j,k)和其水平向上任意一个相邻网格单元说明相邻网格单元处于不同水力联系状态时单元间水平向渗流量的计算方式，用分别表示网格单元(i,j,k)和该相邻网格单元上次迭代求解的潜水位(L)，用/>分别表示本次迭代待求解的潜水位(L)，用Boti,j,k、Botadjoin分别表示底板高程(L)，用C表示网格单元(i,j,k)和该相邻网格单元之间的水力传导系数(L2T-1)：
(1)当时，若/>(图4)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系，本次迭代中渗流量符合传统达西定律，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(2)当时，若/>(图5)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自网格单元(i,j,k)流向该相邻网格单元，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(3)当Boti,j,k≥Botadjoin时，若(图6)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自该相邻网格单元流向网格单元(i,j,k)，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(4)当Boti,j,k≥Botadjoin时，若且/>(图7)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系，单元间无渗流量：
q＝0 (6)
(5)当时，若/>(图8)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系，本次迭代中渗流量符合传统达西定律，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(6)当时，若/>(图9)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自该相邻网格单元流向网格单元(i,j,k)，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(7)当Botadjoin＞Boti,j,k时，若(图10)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自网格单元(i,j,k)流向该相邻网格单元，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(8)当Botadjoin＞Boti,j,k时，若且/>(图11)，则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系，单元间无渗流量：
q＝0 (10)
进一步地，在层方向上，任意两个相邻潜水含水层模拟网格单元之间的渗流量始终符合达西定律，不失一般性，本次迭代中任意一个网格单元(i,j,k)和其沿层方向上的相邻网格单元(i,j,k-1)、(i,j,k+1)之间垂向渗流量分别为：
其中，Ci,j,k-1/2、Ci,j,k+1/2分别表示网格单元(i,j,k)和(i,j,k-1)、(i,j,k+1)之间的水力传导系数(L2T-1)，分别表示网格单元(i,j,k-1)和(i,j,k+1)本次迭代待求解的潜水位(L)。
进一步地，本次迭代中作用于网格单元(i,j,k)之上的源汇项QSi,j,k可由下式计算，
其中，Pi,j,k代表作用于网格单元(i,j,k)之上的与本次迭代待求解潜水位有关的源汇项相关系数(L2T-1)；Qi,j,k代表作用于网格单元(i,j,k)之上的流量源汇项(L3T-1)。
进一步地，本次迭代计算的网格单元(i,j,k)贮水量的变化率为，
其中，SSi,j,k代表网格单元(i,j,k)的贮水率(L-1)，Vi,j,k代表网格单元(i,j,k)的体积(L3)，代表本次迭代计算的网格单元(i,j,k)的潜水位变化率(L1T-1)。
S3、判断本次迭代下各模拟网格单元的模拟潜水位是否满足预设精度阈值，若满足则结束模拟，否则返回S2进行下一次迭代。
进一步地，步骤S3具体为，正常情况下，相邻两次迭代计算的潜水位之差的绝对值(|hm-hm-1|)将逐渐变小，最终近乎相等。当在所有网格单元上本次迭代和上次迭代计算的潜水位相差(|hm-hm-1|)都小于事先规定的精度阈值(hclose)时，认为模拟收敛，即潜水数值仿真矩阵方程求解完成，不再进行迭代，hm即为模拟时段末的模拟潜水位，如图13所示。
本实施例应用实例：
某潜水含水层总面积64km2，含水层一天的降水入渗补给量仅为0.233721m3,含水层底板高程介于4-80米之间，是典型具有大坡度地形条件的干旱区潜水含水层(图12)。潜水自东向西沿底板高程降低方向流出该含水层，要求对其稳定状态下的潜水位进行模拟。行、列方向上网格单元以100m为间距进行剖分，共剖分80行80列，层方向上划分为三层网格单元，共计正方形网格单元19200个(图2)。作为对比，亦尝试使用MODFLOW-2005模拟本实施例的含水层。
本实施例所达到的应用效果：
受大坡度地形条件影响，本实施例水平方向上(行、列方向)相邻网格单元之间很容易出现具有不完全水力联系或完全不具有水力联系的情况，此时MODFLOW-2005不合理地应用传统达西定律计算相邻单元间的水平向渗流，导致模型无法收敛，模型模拟失败；而本发明考虑了水平向渗流修正处理，模型能轻松收敛，本发明模拟的等潜水位线如图13所示。本发明模拟的潜水含水层水量平衡结果如表1所示，相对水量平衡误差可以忽略不计，表明本发明对实施例的模拟是成功的。
表1(单位：m3/d)
以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。

Claims (10)

1.一种大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，包括：
S1、对潜水含水层进行网格化空间离散，获取潜水含水层模拟网格单元；
S2、构建潜水数值仿真矩阵方程，求解本次迭代下各所述模拟网格单元的模拟潜水位；
S3、判断本次迭代下各所述模拟网格单元的模拟潜水位是否满足预设精度阈值，若满足则结束模拟，否则返回S2进行下一次迭代。
2.根据权利要求1所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，对潜水含水层进行网格化空间离散包括：对潜水含水层从行方向、列方向和层方向进行网格化空间离散。
3.根据权利要求1所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，构建所述潜水数值仿真矩阵方程包括：
基于水量平衡关系，构建描述所述模拟网格单元处潜水运动的水量平衡方程；其中，所述水量平衡关系为：在模拟时段内自周围六个相邻网格单元流入模拟网格单元(i,j,k)的潜水量与作用于模拟网格单元(i,j,k)之上的源汇项之和等于模拟网格单元(i,j,k)贮水量的变化；
联立潜水含水层中全部所述模拟网格单元的水量平衡方程，获取所述潜水数值仿真矩阵方程。
4.根据权利要求3所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，所述水量平衡方程为：
其中，qi,j-1/2,k、qi,j+1/2,k、qi-1/2,j,k、qi+1/2,j,k、qi,j,k-1/2、qi,j,k+1/2分别表示网格单元(i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k)、(i,j,k-1)、(i,j,k+1)和(i,j,k)之间的渗流量，渗流量为正值表示潜水流向网格单元(i,j,k)，渗流量为负值表示潜水流出网格单元(i,j,k)；QSi,j,k为作用于网格单元(i,j,k)上的源汇项；表示网格单元(i,j,k)贮水量的变化率，(i,j,k)表示处于第i行、第j列、第k层的模拟网格单元。
5.根据权利要求3所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，所述潜水数值仿真矩阵方程为：
[A]{h}＝{q}
其中，[A]为系数矩阵；{h}为待求解的潜水位向量；{q}为右端项向量，代表水量平衡方程中的常数项和已知项。
6.根据权利要求3所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，求解本次迭代下各所述模拟网格单元的模拟潜水位包括：
在水平方向上，基于预设方程组，对模拟网格单元(i,j,k)及任一相邻网格单元((i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k))之间的水平向渗流量进行计算；其中，(i,j-1,k)、(i,j+1,k)分别为沿行方向上左、右侧相邻的网格单元；(i+1,j,k)、(i-1,j,k)分别为沿列方向上前、后侧相邻的网格单元；
在层方向上，基于预设方程组，对模拟网格单元(i,j,k)及任一相邻的网格单元((i,j,k-1)、(i,j,k+1))之间的垂向渗流量进行计算；其中，其中(i,j,k-1)、(i,j,k+1)分别为沿层方向上上、下侧相邻的网格单元；
将水平方向上各所述模拟网格单元和对应相邻网格单元之间的渗流量、层方向上各所述模拟网格单元和对应相邻网格单元的之间的渗流量、各所述模拟网格单元的源汇项、各所述模拟网格单元贮水量的变化率代入潜水数值仿真矩阵方程求解获取本次迭代下各所述模拟网格单元的模拟潜水位。
7.根据权利要求6所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，在水平方向上，基于预设方程组，对模拟网格单元(i,j,k)及任一相邻的网格单元((i,j-1,k)、(i,j+1,k)、(i-1,j,k)、(i+1,j,k))之间的渗流量进行计算包括：
用分别表示网格单元(i,j,k)和该相邻网格单元上次迭代求解的潜水位，用/>分别表示本次迭代待求解的潜水位，用Boti,j,k、Botadjoin分别表示底板高程，用C表示网格单元(i,j,k)和该相邻网格单元之间的水力传导系数；
(1)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系，本次迭代中渗流量符合传统达西定律，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(2)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自网格单元(i,j,k)流向该相邻网格单元，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(3)当Boti,j,k≥Botadjoin时，若则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自该相邻网格单元流向网格单元(i,j,k)，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(4)当Boti,j,k≥Botadjoin时，若且/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系，单元间无渗流量：
q＝0
(5)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有完全水力联系，本次迭代中渗流量符合传统达西定律，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(6)当时，若/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自该相邻网格单元流向网格单元(i,j,k)，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(7)当Botadjoin＞Boti,j,k时，若则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间具有不完全的水力联系，此时潜水自网格单元(i,j,k)流向该相邻网格单元，但渗流量与/>无关，本次迭代中渗流量需采用修正后的达西定律计算，此时网格单元(i,j,k)和该相邻网格单元之间的渗流量为：
(8)当Botadjoin＞Boti,j,k时，若且/>则认为本次迭代中网格单元(i,j,k)和该相邻网格单元之间完全不具有水力联系，单元间无渗流量：
q＝0。
8.根据权利要求6所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，本次迭代中网格单元(i,j,k)和沿层方向上的相邻网格单元(i,j,k-1)、(i,j,k+1)之间的垂向渗流量分别为：
其中，Ci,j,k-1/2、Ci,j,k+1/2分别表示网格单元(i,j,k)和(i,j,k-1)、(i,j,k+1)之间的水力传导系数，分别表示网格单元(i,j,k-1)和(i,j,k+1)本次迭代待求解的潜水位。
9.根据权利要求6所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，本次迭代作用于所述模拟网格单元之上的源汇项为：
其中，Pi,j,k代表作用于网格单元(i,j,k)之上的与本次迭代待求解潜水位有关的源汇项相关系数；Qi,j,k代表作用于网格单元(i,j,k)之上的流量源汇项；
本次迭代计算的所述模拟网格单元贮水量的变化率为：
其中，SSi,j,k代表网格单元(i,j,k)的贮水率，Vi,j,k代表网格单元(i,j,k)的体积，代表本次迭代计算的网格单元(i,j,k)的潜水位变化率。
10.根据权利要求6所述的大坡度地形条件下的干旱区潜水位模拟方法，其特征在于，判断本次迭代模拟的潜水位是否满足预设精度阈值包括：
将本次迭代模拟的潜水位与上次迭代模拟的潜水位进行比较，当所有所述模拟网格单元上本次迭代模拟的潜水位与上次迭代模拟的潜水位相差都小于预先规定的精度阈值时，认为模拟收敛，即潜水数值仿真矩阵方程求解完成，不再进行迭代。
CN202310512745.5A 2023-05-09 2023-05-09 一种大坡度地形条件下的干旱区潜水位模拟方法 Active CN116522818B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310512745.5A CN116522818B (zh) 2023-05-09 2023-05-09 一种大坡度地形条件下的干旱区潜水位模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310512745.5A CN116522818B (zh) 2023-05-09 2023-05-09 一种大坡度地形条件下的干旱区潜水位模拟方法

Publications (2)

Publication Number Publication Date
CN116522818A true CN116522818A (zh) 2023-08-01
CN116522818B CN116522818B (zh) 2023-12-19

Family

ID=87391884

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310512745.5A Active CN116522818B (zh) 2023-05-09 2023-05-09 一种大坡度地形条件下的干旱区潜水位模拟方法

Country Status (1)

Country Link
CN (1) CN116522818B (zh)

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101908100A (zh) * 2010-07-26 2010-12-08 中国科学院生态环境研究中心 一种地下水环境的建模及数值模拟方法
CN102567634A (zh) * 2011-12-23 2012-07-11 中国水利水电科学研究院 一种基于水循环的地下水数值仿真方法
CN106934093A (zh) * 2017-01-17 2017-07-07 南京大学 模拟三维地下水流运动的三重网格多尺度有限元方法
CN108021780A (zh) * 2018-01-24 2018-05-11 国电南瑞科技股份有限公司 一种基于无规则非结构网格模型的山洪动态仿真方法
CN108763798A (zh) * 2018-06-04 2018-11-06 中国水利水电科学研究院 一种湖泊与地下水非稳定流作用模拟方法
CN110851969A (zh) * 2019-11-01 2020-02-28 中国人民解放军63653部队 水文地质特征与区域地下水循环模拟方法
CN111507026A (zh) * 2019-09-03 2020-08-07 河海大学 一种模拟节点达西渗透流速的双重网格多尺度有限单元法
CN111723505A (zh) * 2020-06-19 2020-09-29 中国水利水电科学研究院 一种流域水质水量监测系统
CN112347678A (zh) * 2020-11-12 2021-02-09 河海大学 一种同时模拟地下水流和达西速度的新型多尺度有限元法
CN112364543A (zh) * 2020-11-17 2021-02-12 中铁第一勘察设计院集团有限公司 基于软塑黄土隧道双排式群井降水模型的渗流模拟方法
US20210181375A1 (en) * 2019-12-12 2021-06-17 China Institute Of Water Resources And Hydropower Research Numerical method for simulating a karez well in association with a groundwater model
CN114547957A (zh) * 2022-02-21 2022-05-27 河海大学 一种结合空间显式水文模型模拟土壤氮淋失的方法
CN115060875A (zh) * 2022-06-17 2022-09-16 大连理工大学 一种基于达西定律的水合物储层生产压力区间的确定方法

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101908100A (zh) * 2010-07-26 2010-12-08 中国科学院生态环境研究中心 一种地下水环境的建模及数值模拟方法
CN102567634A (zh) * 2011-12-23 2012-07-11 中国水利水电科学研究院 一种基于水循环的地下水数值仿真方法
CN106934093A (zh) * 2017-01-17 2017-07-07 南京大学 模拟三维地下水流运动的三重网格多尺度有限元方法
CN108021780A (zh) * 2018-01-24 2018-05-11 国电南瑞科技股份有限公司 一种基于无规则非结构网格模型的山洪动态仿真方法
CN108763798A (zh) * 2018-06-04 2018-11-06 中国水利水电科学研究院 一种湖泊与地下水非稳定流作用模拟方法
CN111507026A (zh) * 2019-09-03 2020-08-07 河海大学 一种模拟节点达西渗透流速的双重网格多尺度有限单元法
CN110851969A (zh) * 2019-11-01 2020-02-28 中国人民解放军63653部队 水文地质特征与区域地下水循环模拟方法
US20210181375A1 (en) * 2019-12-12 2021-06-17 China Institute Of Water Resources And Hydropower Research Numerical method for simulating a karez well in association with a groundwater model
CN111723505A (zh) * 2020-06-19 2020-09-29 中国水利水电科学研究院 一种流域水质水量监测系统
CN112347678A (zh) * 2020-11-12 2021-02-09 河海大学 一种同时模拟地下水流和达西速度的新型多尺度有限元法
CN112364543A (zh) * 2020-11-17 2021-02-12 中铁第一勘察设计院集团有限公司 基于软塑黄土隧道双排式群井降水模型的渗流模拟方法
CN114547957A (zh) * 2022-02-21 2022-05-27 河海大学 一种结合空间显式水文模型模拟土壤氮淋失的方法
CN115060875A (zh) * 2022-06-17 2022-09-16 大连理工大学 一种基于达西定律的水合物储层生产压力区间的确定方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
徐宗学;罗睿;: "PDTank模型及其在三川河流域的应用", 北京师范大学学报(自然科学版), no. 03 *
肖先煊;许模;蔡国军;郭健;吴礼舟;: "基于潜水渗流模型的地下水实际流速", 实验室研究与探索, no. 04 *
赵文智;周宏;刘鹄;: "干旱区包气带土壤水分运移及其对地下水补给研究进展", 地球科学进展, no. 09 *

Also Published As

Publication number Publication date
CN116522818B (zh) 2023-12-19

Similar Documents

Publication Publication Date Title
Mantoglou et al. Management of coastal aquifers based on nonlinear optimization and evolutionary algorithms
Xu et al. Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas
El Yaouti et al. Modelling groundwater flow and advective contaminant transport in the Bou-Areg unconfined aquifer (NE Morocco)
Yihdego et al. Simulation of lake–aquifer interaction at Lake Naivasha, Kenya using a three-dimensional flow model with the high conductivity technique and a DEM with bathymetry
Wu et al. Simulation of soil loss processes based on rainfall runoff and the time factor of governance in the Jialing River Watershed, China
Ahmed et al. Groundwater flow modelling of Yamuna-Krishni interstream, a part of central Ganga Plain Uttar Pradesh
Sivakumar Management policy of water table in dry zone of Sri Lanka to subsidise the pain of non rice crop cultivators for the food productivity improvement.
CN106294282A (zh) 黑油油藏模拟方法及装置
CN110381524A (zh) 基于Bi-LSTM的大场景移动流量在线预测方法、系统及存储介质
CN111090902B (zh) 一种基于地下水模型的坎儿井数值模拟方法
CN105389469A (zh) 一种暴雨洪水管理模型参数自动率定方法
CN110008439A (zh) 基于矩阵分解的降雨数据时空一体化插值算法
Mehl et al. Three-dimensional local grid refinement for block-centered finite-difference groundwater models using iteratively coupled shared nodes: a new method of interpolation and analysis of errors
Saghi-Jadid et al. Result-based management approach for aquifer restoration problems using a combined numerical simulation–parallel evolutionary optimization model
CN116522818B (zh) 一种大坡度地形条件下的干旱区潜水位模拟方法
Kumar Assessing the impact of climate change on groundwater resources
Huntington et al. Integrated hydrologic modeling of Lake Tahoe and Martis Valley mountain block and alluvial systems, Nevada and California
Gerey et al. Groundwater Single‐and Multiobjective Optimization Using Harris Hawks and Multiobjective Billiards‐Inspired Algorithm
Wu et al. The prediction of foundation pit based on genetic back propagation neural network
Jutla Hydrologic modeling of reconstructed watersheds using a system dynamics approach
Bougher et al. Development and validation of the Ground-to-Exosphere Mars GITM Code: Solar cycle and seasonal variations of the upper atmosphere
CN116992783B (zh) 一种地下水大降深疏干下的全有效网格单元流场模拟方法
Ansari et al. Assessment of current and future groundwater stress through varied scenario projections in urban and rural environment in parts of Meerut district, Uttar Pradesh in Ganges sub-basin
Carroll et al. An unconfined groundwater model of the Death Valley Regional Flow System and a comparison to its confined predecessor
Gallardo et al. Groundwater supply under land subsidence constrains in the Nobi Plain

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
GR01 Patent grant
GR01 Patent grant