1、水文地质数值计算水文地质数值计算 主要参考书目：1 地下水流动问题数值方法，陈崇希等，1990.6，中国地质大学出版社。2 地下水流数值模拟，1989.9,李俊亭，地质出版社 3 地下水动力学（第二版），薛禹群，1997.9，地质出版社。4 地下水流的数学模型和数值法，孙讷正，1981.5，地质出版社。教学中涉及软件：1 Processing Modflow(ETH)2 Surfer(V.6或更高),Golden Software Inc.3 Auto CAD 2000(Autodesk Inc.)4 Mcrosoft Excel达西定律微分表达式1 渗流理论基础知识（1）1.1达西公式及其微
2、分表达式21,xyzxxxyyyzzzHHdHvQKAvKJvKuLdsndHdHdHvKvKvKdxdydzdHdHdHvKvKvKdxdydz 各向同性介质各向异性介质，主方向与所选座标方向平行 xxxxyxzyyxyyyzzzxzyzzdHdHdHvKKKdxdydzdHdHdHvKKKdxdydzdHdHdHvKKKdxdydz 各向异性介质，主方向与所选座标方向平行 1 渗流理论基础知识（2）1.2 地下水流动微分方程地下水流动微分方程-基本微分方程基本微分方程各 向 同 性 介 质)物 理 含 义？-d i v(v)量 纲？含 义？（又 可 写 为）量 纲？含 义？估 算 公 式
3、为“渗问 题杂”流中 最 复的 方 程S()()()()()()(SSSSHHHHKKKWSxxyyzztHHHKKKxxyyzzWSSSgnS=-=-每 降 低 一 个 单 位 压 强，单 位 体 积 水 所 膨 胀 增 加 的 体 积每 降 低 一 个 单 位 压 强，单 位 体 积 的 地 层 压 缩“挤”出 水 的 体 积单 位 体 积 的 含 水 层 中（包 括 水 和 岩 石“骨 架”），每 降 低 一 个 单 位 水 头，所 释 放 出 水 的 体介 质，当选标与质渗平 行 时积对 于 各 向 异 性所座方 向介主方 向)11 -()(bbx xd Vd VVd pVd pSHK
4、Kxxy)()yyzzSHHHKWSyzzt思考问题：“地下水流动微分方程地下水流动微分方程”代数方程的解为未知数代数方程的解为未知数(x1,x2),x1,x2),微分方程的解是什么？微分方程的解是什么？未知函数为什么建立微分方程为什么建立微分方程？可求解出水头的空间分布，进而可解出水压强、流速、流量等为什么建立关于水位为什么建立关于水位“H”H”的微分方程？的微分方程？比压强P、流速v等物理量的方程简单 ()-()SHHHg ra d HijkxyxvKg ra d HHd iv vWSt达 西 定 律(各 向 同 性 介 质)梯 度 概 念？微 分 方 程（散 度 概 念？）渗 流 问 题
5、 的 场 论 表 示 法：1 渗流理论基础知识（2）1.2 地下水流动微分方程地下水流动微分方程-承压水流的微分方程式承压水流的微分方程式 ()()()1 SSSS对 基 本 微 分 方 程 沿 垂 向 积 分，上 下 界 取 含 水 层 顶 底 板，为 含 水 层 内 垂 向 平 均 水 头（流 动 方 向 主 要 为 水 平 分 量）注 意：此 处的 含 义？量 纲？、的 含 义？量 纲？估 算 公 式 为=-=HHHHTTWSxxyytWSMSSSSSgndVVdp 最最 简简 单单 的的 单单 层层 承承 压压 含含 水水 层层 渗渗 流流 方方 程程1 -S-每 降 低 一 个 单
6、位 压 强，单 位 体 积 水 所 膨 胀 增 加 的 体 积每 降 低 一 个 单 位 压 强，单 位 体 积 的 地 层 压 缩“挤”出 水 的 体 积单 位 体 积 的 含 水 层 中（包 括 水 和 岩 石“骨 架”），每 降 低 一 个 单 位 水 头，所 释 放 出 水 的 体 积bbdVVdpS1 渗流理论基础知识（2）1.2 地下水流动微分方程地下水流动微分方程 潜水流的微分方程式潜水流的微分方程式()()yyHHHHK HBK HBWSxxytWS 对基本微分方程沿垂向积分，上界取潜水面，下界取潜水含水层底板，为含水层内垂向平均水头注意：此处的含义？量纲？的含义？量纲？(或)
7、的含义？量纲？（引入裘布依假定）最最简简单单的的单单层层潜潜水水含含水水层层“近近似似”微微分分方方程程 非非线线性性方方程程！1 渗流理论基础知识（2）1.2 地下水流动微分方程地下水流动微分方程 多层含水层越流系统的微分方程式多层含水层越流系统的微分方程式准三维流准三维流 之一（第二类越流系统）之一（第二类越流系统）1111112222()()yHHHK HBK HBWSxxytHHHTTWSxxyt 以两层为例，上潜水、下承压地下水流系统。含水层内主要为水平流动分量，弱透水层内主要为垂直流动分量潜水：承压水：忽略弱透水层的弹性与压缩释放水（第二类越流系统）2 21 11 12 2K K(
8、H H-H H)mmK K(H H-H H)mm尽尽管管用用两两个个（或或多多个个）二二维维渗渗流流方方程程联联立立，但但还还是是比比三三维维问问题题简简单单！1 渗流理论基础知识（2）1.2 地下水流动微分方程地下水流动微分方程 多层含水层越流系统的微分方程式多层含水层越流系统的微分方程式准三维流准三维流 之二（第三类越流系统）之二（第三类越流系统）11211111122333()()yZSHHHHKHBKHBWKSxxyztHHKSzztHHTTWKxxy 潜水底板以两层为例，上潜水、下承压,含水层内主要为水平流动分量，弱透水层内主要为垂直流动分量潜水：弱透水层承压水：考虑弱透水层的弹性释
9、水（第三类越流系统）32ZHHSzt承压顶板三，尽尽管管用用个个（或或多多个个）二二维维渗渗流流方方程程联联立立仍仍然然比比三三维维问问题题简简单单！思考问题：“地下水二维与准三维流微分方程地下水二维与准三维流微分方程”真三维流地下水方程什么样？真三维流地下水方程什么样？基本微分方程在二维、准三维流含水层系统中，含水层内的水头沿垂直方向变化否？在二维、准三维流含水层系统中，含水层内的水头沿垂直方向变化否？不变为什么准三维流比真三维流的方程还复杂、麻烦？为什么准三维流比真三维流的方程还复杂、麻烦？其实是简单了，少了一个空间变量第二类越流系统中，弱透水层内的水头是怎样分布的？第二类越流系统中，弱透
10、水层内的水头是怎样分布的？连接上下含水层水头的直线一般情况下，弱透水层的储水率大还是含水层的储水率大？一般情况下，弱透水层的储水率大还是含水层的储水率大？弱透水层的储水率要比含水层大，地面沉降主要发生在弱透水层 边界条件边界条件 地下水系统或模型边界的水头、流量。（1）第一类边界，A 已知水头边界：hA=h1(x,t)（2）第二类边界，B已知流量边界：隔水边界（C）流量为零的二类边界。初始条件初始条件非稳定流问题非稳定流问题 t=0时地下水系统或模型的水头分布。QABCh1(,)BhkMq x y tn定解问题：“地下水流泛定方程”定解条件泛定方程：不带边界条件与初始条件的微分方程定解条件：又
11、分为边界条件、初始条件等。与泛定方程联立后，可使其具有唯一解。稳定流问题没有初始条件！边界条件：边界条件：对所求解微分方程的未知函数未知函数,给出边界上的有关已知信息已知信息，常用有第一类边界条件（已知已知水位边界）、第二类边界条件 （已知已知流量边界已知已知未知函数边界法向导数）等。对于稳定流问题，如全部为已知流量边界（第二类边界）条 件是不适定的（无唯一解）。即至少有一段（片）为一类边 界条件 初始条件：初始条件：微分方程求解的是某初始时刻初始时刻以后的变化过程。初始时刻是选定的，因问题与目的而异。初始条件要求给出：初始时刻的未知函数的数值或状态。1 渗流理论基础知识（1）1.3 地下水流
12、定解问题边界形状边界形状：三维问题边界：封闭曲面（可多个组合）二维问题边界：封闭曲线（可多个组合）一维问题边界：两个端点对于稳定流问题：至少有一段对于稳定流问题：至少有一段(点点)是第一类边界条件，否则是无解是第一类边界条件，否则是无解的。的。为什么仅归纳出三种边界类型？为什么仅归纳出三种边界类型？还有其它类型的边界条件吗？还有其它类型的边界条件吗？，已知函数或数值i ii ii ii ip pp pp pp pp pp pp pp p-2 2p pp p开开采采强强度度 函函数数 表表示示法法（脉脉冲冲函函数数）Q Q (x x-x x,y y-y y)函函数数特特征征1 1 (x x-x
13、x,y y-y y)d dx xd dy y=1 1 ,x x=x x,y y=y y2 2 (x x-x x,y y-y y)=0 0 ,x x x x,y y y y3 3 (x x-x x,y y-y y)量量纲纲 L L 矩形网格有限差剖分图矩形网格有限差剖分图通过结点编号间接表示结点的位置：单编号法：简单，但位置不明确 双编号法：复杂，但位置隐含其中（推荐）矩形网格有限差节点编号方法矩形网格有限差节点编号方法i,j,m,ki,j,m,kmmi,j,ki,j,ki,j,ki,j,k三三维维问问题题H H 含含义义准三准三维维问问题题H H 含含义义二二维维问问题题H H 含含义义,1,
14、1/2,1,1/2,1,1/2,1,1/2,1,1/2()()()()ij kij kijkijij kij kijij kij kijijkij kijijkij kij kHHSxytxTHHTHHyyTHHTHHWxyx矩形网格显式差分格式矩形网格显式差分格式（五点格式），可用水均衡法导出。（五点格式），可用水均衡法导出。特点：每个结点的方程中仅仅有一个未知数，不用联立即可求解。-2-2222212()()()1()4()4()1()2()2StM inxyTxyxSTttM inrTxSxSTttM inrTxS 当当时时，或或 二二 维维 流流 情情 况况 或或 一一 维维 流流 情
15、情 况况,1,1/2,1,2i jijiji jijT TTTT稳定与收敛条件稳定与收敛条件收敛性：收敛性：在计算中没有舍入误差的情况下在计算中没有舍入误差的情况下，当x、y、t、趋于零时，差分方程的解hi,j收敛于微分方程的解H(x,y,t)，即截断误差截断误差趋近 于零（用差商代替微商时，引入了截断误差），称该差分方程 是收敛的。稳定性稳定性：在有一定舍入误差或概化误差情况下，当x、y、t、趋于零，差分 方程的“实际”解Hi,j收敛于差分方程的“精确解”hi,j称该差分方程（格式）是稳定的。显式差分格式是有显式差分格式是有条件稳定与收敛的条件稳定与收敛的一维流情况二维流情况时，当 2)(4
16、)(yx)()(212222TSxMintTSxMintTSyxMint微分方程解微分方程解=差分方程解差分方程解=（有舍入与采集误差的）代数方程解（有舍入与采集误差的）代数方程解收敛性收敛性稳定性稳定性（对于非稳定递推问题）（对于非稳定递推问题）一维问题计算例22,01000,0(,0)1(0,)0,(1000,)0hhTSxtxth xhtht定解问题描述：221000/,100,0.01,0.01 TmdxmStdTtrx S 符符合合一一维维流流计计算算稳稳定定与与收收敛敛条条件件1 11 1 1 10 02 2右侧右侧数值计算结果数值计算结果221000/,100,0.01,0.04
17、8 TmdxmStdTtrx S 符符合合一一维维流流计计算算稳稳定定与与收收敛敛条条件件1 10 0.4 48 8 2 2221000/,100,0.01,0.1 TmdxmStdTtrx S不符符合合一一维维流流计计算算稳稳定定与与收收敛敛条条件件1 11 1 2 2显式计算格式操心两件事：显式计算格式操心两件事：1 x、y、t 不能太大，否则不收敛。不能太大，否则不收敛。2 t与与x、y之间的关系要满足一定的约束，否则不稳定之间的关系要满足一定的约束，否则不稳定 若不满足下面的约束式，若不满足下面的约束式，DtDt再小也不收敛！再小也不收敛！收敛与稳定条件：收敛与稳定条件：-2-2222
18、212()()()1()4()4()1()2()2StMinxyTxyxSTttMinrTxSxSTttMinrTxS 当当时时，或或 二二 维维 流流 情情 况况 或或 一一 维维 流流 情情 况况,1,11,1,111,1,1,221,1,1,11,1,1,1,221,2,1,1,2()()()()2 i j ki j ki jkij ki j kij ki j kijiji jki j ki jki j ki ji ji j ki jijijiHHSxytxTHHTHHyyTHHTHHxWxyTTTT1,1,2 ()调和平均，符合非均质平均渗透系数计算方法）含义？jiji j kTW二维
19、流简隐式差分格式二维流简隐式差分格式特点特点：方程中五个未知数，必须联立才能求解！方程中五个未知数，必须联立才能求解！简隐式差分格式无条件稳定与收敛。简隐式差分格式无条件稳定与收敛。对时间步长不受稳定条件约束限制对时间步长不受稳定条件约束限制,1,1/2,1,1,11/2,1,1,1,1/2,1,1,1,1/2,1,1,11/2,1,1/2()()()()()ij kij kijkijij kij kijij kij kijijkij kijijkij kijij kij kiHHSxytxTHHTHHyyTHHTHHxTHHT,1,1/2,1,1/2,1,1/2()()()*0.5 jij
20、kij kijijkij kijijkij kij kxHHyyTHHTHHxWxy二维流中心（对称）差分格式二维流中心（对称）差分格式特点特点：方程中五个未知数，必须联立才能求解！方程中五个未知数，必须联立才能求解！无条件稳定与收敛（对时间步长不受稳定条件约束限制）。无条件稳定与收敛（对时间步长不受稳定条件约束限制）。稳定性没有简隐式好，收敛性超过简隐式。稳定性没有简隐式好，收敛性超过简隐式。实际问题多取简隐式。实际问题多取简隐式。1 那些结（格）点列出差分方程式那些结（格）点列出差分方程式 内结（格）点=列出差分式 流量边界（第二类）结（格）点=列出差分式 水位边界（第一类）结（格）点=不
21、列出差分式 外结（格）点=不列出差分式2 所有内结点的方程式都是相同的（给编写程序带来方便）3 当导水系数用调和平均值时（命：计算区域之外的导水系数为零），隔水边界（零流量边界）的差分式与内结点相同。4 非零流量边界“等价于”隔水边界+“注水井”，注水井流量的大小 等于边界流入量。,1,11,1,111,1,1,221,1,1,11,1,1,1,22 ()()()()i j m ki j m kS i jkij m ki j m kij m ki j m kij mij mi jm ki j m ki jm ki j m ki jmi jmHHSx y zty zkHHkHHyx zkHHkH
22、Hx1,1,1,11,1,1,1,221,2()()x y z i j mki j m ki j mki j m ki j mi j mi j m kx ykHHkHHzW三维问题的有三维问题的有“简隐式简隐式”限差分方程限差分方程 相对弱透水的相对弱透水的“夹层夹层”也作为计算层处理。也作为计算层处理。即不仅计算含水层水位，同时计算弱透水层中的水位。即不仅计算含水层水位，同时计算弱透水层中的水位。,1,11,1,111,1,1,221,1,1,11,1,1,1,2211,1,2,()()()()(mmi j ki j kmi jkmmmmmmij ki j kij ki j kijijmmm
23、mmmi jki j ki jki j ki ji jmmmi j ki j ki jHHSx ytxTHHTHHyyTHHTHHxHHk111,1,12,1/2)()1,2,3.mmmi j ki j ki ji j kHHkx ymmWx ym准三维问题的准三维问题的“简隐式简隐式”有限差分方程有限差分方程 仅仅计算各个含水层中的水位，其他仅仅计算各个含水层中的水位，其他“夹层夹层”不计算不计算 忽略夹层中的弹（塑）性储存，将其近似地概化化到相邻的含水层中去。忽略夹层中的弹（塑）性储存，将其近似地概化化到相邻的含水层中去。稳定流数学模型与非稳定流之区别：稳定流数学模型与非稳定流之区别：1
24、稳定流模型是非稳定流模型的特例稳定流模型是非稳定流模型的特例2 稳定流模型中不要求初始条件稳定流模型中不要求初始条件3 稳定流模型中，边界条件不能全部是流量边界（二类边界），至少稳定流模型中，边界条件不能全部是流量边界（二类边界），至少 有一点（段）是一类或三类边界。有一点（段）是一类或三类边界。4 方程中若方程中若 并不一定是稳定流模型，边界条件中有非稳定因素仍然并不一定是稳定流模型，边界条件中有非稳定因素仍然 是非稳定流模型。是非稳定流模型。0Ht0Htv其它剖分类型（二维）其它剖分类型（二维）三角形网络三角形网络四边形网络四边形网络v离散离散网格的生成（三维）网格的生成（三维）v离散离散
25、网格的生成（三维）网格的生成（三维）立体剖分立体剖分 v离散离散网格的生成（三维）网格的生成（三维）立体剖分立体剖分 拱坝拱坝前处理：参数分布前处理：参数分布1234非非均质问题的近似描述：均质问题的近似描述：分片常数法分片常数法前处理：参数分布前处理：参数分布非非均质问题的近似描述：均质问题的近似描述：等值线插值法等值线插值法 点域点域“水均衡水均衡”差分格式差分格式 非对称网格有限差分方程非对称网格有限差分方程p12345iDiLipipipipipiDLhhTQ)()(pipipipipipipiphhChhLDTQQpeppSthSQthhhhSCnpnpenpnipp1)(数值模拟：
26、特殊问题的处理 格点水位与开采井水位流量和水头数值模拟：特殊问题的处理 非线性问题解法线性化tHHxHHKx)()(潜水非饱和流可压缩含水层吸附性溶质运移t=t+1H0,K0,0k=0k=k+1Hk+1=#(Hk,Kk,k)Kk+1=K(Hk+1)k+1=(Hk+1)迭代误差数值模拟：特殊问题的处理 潜水自由面(浸润面)干湿迭代差分法干湿迭代差分法移动网格有限元法移动网格有限元法固定网格有限元法固定网格有限元法修正渗透矩阵法初流量法分裂法数值模拟：特殊问题的处理 数值弥散和振荡ixCVxCDtC22VNaClC0C0 xCxxCCVxCCCDtCCninininininini221111211
27、1111Pelect数：Vx/D数值弥散数值弥散0数值振荡数值振荡数值模拟：特殊问题的处理 数值弥散和振荡0,/)(0,/)(11VxCCVxCCxCiiii迎风格式，上游差分格式迎风格式，上游差分格式上游加权差分法上游加权有限元法数值振荡减小，但数值弥散增加移动坐标系统移动网格特征有限元法：质点跟踪技术随机游动方法：对流+随机质点数值模拟：常用软件 地下水模拟 美国地质调查局(UGS):MODFLOW,MT3D Processing ModFlow(PM),Visual ModFlow,Groundwater Modeling System(GMS)加拿大阿波罗集团(APOLLO)：FEFLOW
