1、地下水数值模拟原理及地下水数值模拟原理及建模方法和步骤建模方法和步骤 内容提要内容提要绪论绪论一、地下水流有限差分法原理一、地下水流有限差分法原理二、地下水数值模型建模步骤二、地下水数值模型建模步骤三、建模所需要的基本资料三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流五、地下水系统热对流-扩散有限差分法原扩散有限差分法原理理绪论绪论地下水数值方法在水文地质学中的位置地下水数值方法在水文地质学中的位置o 已学课程已学课程n 水文地质学基础水文地质学基础n 地下水动力学地下水动力学n 水文地球化学（环境水文地质学）水文地球化学
2、（环境水文地质学）n 同位素同位素n 地下水溶质运移：示踪剂或污染物在含水层中地下水溶质运移：示踪剂或污染物在含水层中的迁移机理及数学模型和求解方法的迁移机理及数学模型和求解方法o 地下水数值模拟地下水数值模拟绪论绪论地下水数值方法在水文地质学中的位置地下水数值方法在水文地质学中的位置地下水动力学主要内容地下水动力学主要内容o连续性原理、达西定律、水均衡原理连续性原理、达西定律、水均衡原理o地下水流基本方程地下水流基本方程o几类特殊水文地质问题数学模型及解析解几类特殊水文地质问题数学模型及解析解n地下水向沟渠河中的流动地下水向沟渠河中的流动n园岛模型园岛模型n泰斯模型泰斯模型n有越流的不稳定井
3、流（有越流的不稳定井流（Hantush and Jacob）n无越流的潜水含水层不稳定井流（无越流的潜水含水层不稳定井流（Neuman）o抽水试验及反求参数抽水试验及反求参数o数值方法可以作为地下水动力学课程内容之一、或补充、或延续数值方法可以作为地下水动力学课程内容之一、或补充、或延续绪论绪论数值解与解析解数值解与解析解o 地下水动力学中所得到的解是解析解地下水动力学中所得到的解是解析解o 解析解的特点解析解的特点n 可以利用解析解公式计算出研究区任意空间点可以利用解析解公式计算出研究区任意空间点和时间点处的水头值和时间点处的水头值n 对实际模型的概化对实际模型的概化n 适用于一些特殊简单问
4、题适用于一些特殊简单问题n 公式复杂，需要借助计算机求解公式复杂，需要借助计算机求解o 因此，其应用受到很大限制因此，其应用受到很大限制绪论绪论数值解与解析解数值解与解析解o 数值方法是地下水动力学的完善和补充或延续数值方法是地下水动力学的完善和补充或延续o 数值解的特点：数值解的特点：n只是求出研究区某些空间点和某些时间点处的水头值只是求出研究区某些空间点和某些时间点处的水头值n适用于所有的问题适用于所有的问题n具备水文地质基础和线性代数知识具备水文地质基础和线性代数知识n已有数值模拟专门软件（或自己编程）已有数值模拟专门软件（或自己编程）n需要有高性能计算机需要有高性能计算机n对实际问题的
5、刻画比较精确对实际问题的刻画比较精确o 因此，其应用非常广泛因此，其应用非常广泛数值方法很多，但是最简单实用的是有限差分法：数值方法很多，但是最简单实用的是有限差分法：绪论绪论数值方法数值方法 有限差分法有限差分法有限单元法有限单元法积分有限差分法积分有限差分法半解析半数值法半解析半数值法边界元法边界元法有限体积法有限体积法只讲有限差分法只讲有限差分法一、有限差分法的基本原理一、有限差分法的基本原理有两种方法建立差分方程有两种方法建立差分方程 方法一方法一以地下水流基本微分方程及以地下水流基本微分方程及其定解条件为基础，其定解条件为基础，在渗流区剖分在渗流区剖分基础上，用差商代替微商，将地下水
6、基础上，用差商代替微商，将地下水流微分方程的求解转化为差分方程流微分方程的求解转化为差分方程（代数方程）求解。（代数方程）求解。适用于二维矩形适用于二维矩形网格剖分、三维长方体网格剖分网格剖分、三维长方体网格剖分。方法二在渗流区剖分的基础上，直方法二在渗流区剖分的基础上，直接由达西定律和水均衡原理，建立各接由达西定律和水均衡原理，建立各个均衡区的水均衡方程，从而得到差个均衡区的水均衡方程，从而得到差分方程分方程。适用于矩形网格、三角形网适用于矩形网格、三角形网格。格。矩形网格矩形网格多边形网格多边形网格1、网格划分的基本类型、网格划分的基本类型o（1）先划格线，格点位）先划格线，格点位于网格中
7、心于网格中心均衡网格均衡网格节点网格节点网格o（2）先规定格点位置，）先规定格点位置，再垂直平分两相邻结点的连再垂直平分两相邻结点的连线作格线，形成的网格即为线作格线，形成的网格即为水均衡区水均衡区方法一：差商代替微商方法一：差商代替微商MODFLOW网格系统网格系统方法一：差商代替微商方法一：差商代替微商x000()()()f xxf xfxx 导导数数的的有有限限差差商商近近似似导数的定义导数的定义 当当非常小的时候，有非常小的时候，有 上式右端项即为上式右端项即为f(x)(x)在在x0 0处的差商。处的差商。这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导
8、出差商及其误差。xxfxxfxfx)()(lim)(0000方法一：方法一：差商代替微商差商代替微商(2)有限差分方程建立有限差分方程建立2000()()()()()2!ff xxf xfxxx 2000()()()()()2!ff xxf xfxxx 000()()()()f xxf xfxOxx 已知泰勒公式 由A得：AB 由B 得：000()()()()f xf xxfxOxx 称称 为为f(x)在在x0处的处的一阶前向差商，一阶前向差商，为为截断误差截断误差。xxfxxf)()(00)(xO 称称 为为f(x)在在x0处的处的一阶后向差商，一阶后向差商，为为截断误差截断误差。xxxfx
9、f)()(00)(xO 方法一方法一200002()2()()()()()f xxf xf xxfxOxx 由A-B可以得：由A+B可以得：2000()()()()2f xxf xxfxOxx (4)23400000()()()()()()()()()2!3!4!fxfxff xxf xfxxxxx (4)23400000()()()()()()()()()2!3!4!fxfxff xxf xfxxxxx AB称称 为为f(x)在在x0处的处的一阶中心差商，一阶中心差商，为为截断误差截断误差。xxxfxxf2)()(002)(xO 称称 为为f(x)在在x0处的处的二阶二阶中心中心差商，差商，
10、为为截断误差截断误差。2000)()()(2)(xxxfxfxxf2)(xO 方法一方法一o 对于偏导数（偏微商），类似可以得到相应的差商：对于偏导数（偏微商），类似可以得到相应的差商：ttxHttxHttxH),(),(),(000000 xtxHtxxHxtxH),(),(),(00000020000002002)(),(),(2),(),(xtxxHtxHtxxHxtxH方法一方法一(2)有限差分方程建立有限差分方程建立(续续)一维控制方程差分格式一维控制方程差分格式ttxhxtxhT),(),(22ttxHttxHxtxxHtxHtxxH),(),()(),(),(2),(000020
11、00000ttxHttxHxttxxHttxHttxxH),(),()(),(),(2),(00002000000显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一1,3,2 )(21211 nxitHHxHHHninininini控制方程控制方程1,3,2 )(21211111 nxitHHxHHHninininini网格剖分网格剖分nx个个二维控制方程差分格式二维控制方程差分格式thyhTxhT2222tHHyHHHxHHHnjinjinjinjinjinjinjinji,1,21,1,2,1,1)(2)(2显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一控制方程控制
12、方程网格剖分网格剖分nx个个tHHyHHHxHHHnjinjinjinjinjinjinjinji,1,211,1,11,21,11,1,1)(2)(2达西定律达西定律：水均衡原理水均衡原理：对某一研究对象，流入流入 流出流出体系内质量（或水量）变化量体系内质量（或水量）变化量研究对象可以是大区域的，也可以是微分单元体大区域的水均衡计算经常用于区域的水资源评价dxdHKq单位时间通过单位面积过水断面的水单位时间通过单位面积过水断面的水量与断面处的水力梯度成正比。量与断面处的水力梯度成正比。方法二：直接由达西定律和水均衡原理建立方法二：直接由达西定律和水均衡原理建立差分方程差分方程 表示：当水头
13、下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）e2)(2)(111,1,2/1111,1,2/1jjiijijijijjiijijijiyyxxHHKMyyxxHHKMjijiiijjjijijiiijjjijijiAxxyyHHKMxxyyHHKM,111,1,2/1,111,1,2/1,2)(2)(+H 可取可取tn+1或或 tn时时刻的值刻的值节点（节点（i,j）的均衡区）的均衡区Aij:节点节点(i,j)均衡区的面积均衡区的面积jinnnjinjijeiAttHH,1,1,(i,j)(i+1,j)(i
14、-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例以二维承压水流为例 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）eH取tn时刻的值时刻的值，有：2)(2)(111,1,2/1111,1,2/1jjiinjinjijijjiinjinjijiyyxxHHKMyyxxHHKMjijiiijjnjinjijiiijjnjinjijiAxxyyHHKMxxyyHHKM,111,1,2/1,111,1,2/1,2)(2)(+节点（节点（i,j
15、）均）均衡区衡区Aij:节点节点(i,j)均衡区的面积均衡区的面积jinnnjinjijeiAttHH,1,1,(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续1）jinnjeijiAttC,1,2)(,2)(111,2/1,111,2/1,jjiijijijjiijijiyyxxKMEyyxxKMW2)(,2)(1112/1,1112/1,iijjjijiiijjjijixxyyKMNxxyyKMS(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xyjinj
16、ijijijinjijijijijinjijinjijinjijinjijinjiCHCAHNSEWHNHSHEHWH,1,1,1,1,1,/)(则得到显式格式：则得到显式格式：不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续2）2)(2)(1111,1,1,2/11111,1,1,2/1jjiinjinjijijjiinjinjijiyyxxHHKMyyxxHHKMjijiiijjnjinjijiiijjnjinjijiAxxyyHHKMxxyyHHKM,1111,11,2/1,1111,11,2/1,2)(2)(+如果如果 H 取取tn+1时
17、刻的值时刻的值节点（节点（i,j）的均衡区）的均衡区jinnnjinjijeiAttHH,1,1,(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续3）jinnjeijiAttC,1,2)(,2)(111,2/1,111,2/1,jjiijijijjiijijiyyxxKMEyyxxKMW2)(,2)(1112/1,1112/1,iijjjijiiijjjijixxyyKMNxxyyKMS记：记：jijinjijinjijijijijijinjijinjijinjijinjiji
18、AHCHCNSEWHNHSHEHW,1,11,11,1,1,1,1,)(得到隐式格式：得到隐式格式：(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续4）表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）e11,jjiiyyyxxxyxtHHnjinjie,1,当网格是等距时，即yxHHKMyxHHKMjijijiji,1,1yxxyHHKMxyHHKMjijijijiji,1,1,+H
19、 可取可取tn+1或或 tn时时刻的值刻的值节点（节点（i,j）的均衡区）的均衡区此时，有水均衡方程：(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续5）(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy节点（节点（i,j）的均衡区）的均衡区两边除以 ，得到 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）eyxtHHnjinjie,1,21,1,121,1,1)
20、()(xHHKMxHHKMnjinjinjinjijinjinjinjinjiyHHKMyHHKM,21.11.21.11.)()(H 取取tn+1时刻的值，得时刻的值，得到隐式格式到隐式格式不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续6）表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）etyxe/njinjiHH,1,两边除以21,1,121,1,1)()(xHHtKMxHHtKMnjinjienjinjieejinjinjienjinjiety
21、HHtKMyHHtKM/)()(,21.11.21.11.+H 可取可取tn+1时刻的值，时刻的值，得到隐式格式得到隐式格式节点（节点（i,j）的均衡区）的均衡区(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续7）(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）e如果（i,j)是左边界节点，有2)(2111,1,
22、2/111,jjiijijijijjjiyyxxHHKMyyqjijiiijjjijijiiijjjijijiAxxyyHHKMxxyyHHKM,11,1,2/1,11,1,2/1,2)(2)(+H 可取可取tn+1或或 tn时时刻的值刻的值节点（节点（i,j）的均衡区）的均衡区qij:节点节点(i,j)处的边界单宽流量。处的边界单宽流量。jinnnjinjijeiAttHH,1,1,不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续8）三维有限差分格式（三维有限差分格式（显示显示）222111111,kkjjiikjizzyyxxVkjinkji
23、nkjikjsikjinkjikjikjikjikjikknkjinkjikjzikknkjinkjikjzikjikjikjikjijjnkjinkjikjxijjnkjinkjikjyikjikjikjikjiiinkjinkjikjxiiinkjinkjikjxiVHHVzzxxzzHHKzzHHKzzxxyyHHKyyHHKzzyyxxHHKxxHHK,1,1,1,1,11,1,2/1,1,1,2/1,1,1,1,11,1,2/1,1,1,2/1,1,1,1,1,1,1,2/11,1,2/1222222三维有限差分格式（三维有限差分格式（隐式隐式）222111111,kkjjiikji
24、zzyyxxVkjinkjinkjikjsikjinkjikjikjikjikjikknkjinkjikjzikknkjinkjikjzikjikjikjikjijjnkjinkjikjxijjnkjinkjikjyikjikjikjikjiiinkjinkjikjxiiinkjinkjikjxiVHHVzzxxzzHHKzzHHKzzxxyyHHKyyHHKzzyyxxHHKxxHHK,1,1,1,1,11,11,2/1,1,11,2/1,1,1,1,11,1,1,2/1,1,1,1,2/1,1,1,1,1,1,1,1,2/11,1,1,2/1222222 首先将渗流区划分成若干个辅助小三角
25、形（图36中虚线所示）。划分时应注意：三角形的任一内角不得大于90，三条边的长度尽可能接近；三角形的顶点不能落在另外某个三角形的边上；相邻三角形变化不太大；应考虑水文地质条件，尽量使抽水井、观测孔（特别是用于拟合的观测孔）位于三角形的顶点处，河流、断层等位于三角形边上。1渗流区的剖分渗流区的剖分任意多边形网格系统任意多边形网格系统三角形的顶点称为格点，三角形的顶点称为格点，分为内格点和边界格点分为内格点和边界格点（包括第一类边界格点（包括第一类边界格点和第二类边界格点）。和第二类边界格点）。2、任意多边形网格有限差分法、任意多边形网格有限差分法一、有限差分法的基本原理一、有限差分法的基本原理
26、2多边形均衡区网格的形成方法是：对某个格点i，分别作出格点i与其相邻格点的连线的垂直平分线，这些平分线所围成的多边形区域即为格点 i 的多边形均衡网格子区对每个格点都这样做，则形成多边形均衡网格系统（图36中实线所示）ijkpqacbxyo 典型多边形子区典型多边形子区2、任意多边形网格有限差分法、任意多边形网格有限差分法 ijkpqacbxyo 典型多边形子区典型多边形子区innninieiiieninkikninjijAttHHAqbk iHHTpbj iHHT111111i=1,2,N隐式格式隐式格式1111nnnnjikieijikhhhhQTpbTbqijik1111nnnnjiki
27、ijikehhhhQTpb Tbqijik侧ViiiQA水均衡方程为：水均衡方程为：2、任意多边形网格有限差分法、任意多边形网格有限差分法 ijkpqacbxyo 典型多边形子区典型多边形子区innninieiiieninkikninjijAttHHAqbk iHHTpbj iHHT11i=1,2,N显式格式显式格式2、任意多边形网格有限差分法、任意多边形网格有限差分法o（1 1）区域剖分）区域剖分 首先，在平面上将研究区剖分成若干个辅助小三角形，首先，在平面上将研究区剖分成若干个辅助小三角形，然后以这些三角形为底面，将区域剖分成一些垂直的柱形体，然后以这些三角形为底面，将区域剖分成一些垂直的
28、柱形体，再根据地层岩性及构造特征，用一些平面将柱形体剖分成若再根据地层岩性及构造特征，用一些平面将柱形体剖分成若干短柱体（图干短柱体（图3.13.1）。）。关于平面三角形的剖分依据下述原则：关于平面三角形的剖分依据下述原则：1 1）三角形的任）三角形的任一内角不得大于一内角不得大于9090，三条边的长度尽可能接近；，三条边的长度尽可能接近；2 2）三角形）三角形的顶点不能落在某个三角形的边上；的顶点不能落在某个三角形的边上；3 3）相邻三角形变化不太）相邻三角形变化不太大。大。3 3、多棱柱体三维流有限差分法、多棱柱体三维流有限差分法o（2 2）均衡区网格的形成）均衡区网格的形成 三角形的顶点
29、称为格点，分为内格三角形的顶点称为格点，分为内格点和边界格点（包括第一类边界格点和点和边界格点（包括第一类边界格点和第二类边界格点）。第二类边界格点）。各结点的均衡区形成方法是：平面各结点的均衡区形成方法是：平面上，以格点为中心，与之相连的边的垂上，以格点为中心，与之相连的边的垂直平分线段组成一个多边形；垂向上，直平分线段组成一个多边形；垂向上，以该层与上下两层的两个中间层的平面以该层与上下两层的两个中间层的平面分别为上下底面，这样得到一个以平面分别为上下底面，这样得到一个以平面上为多边形的立体柱体，就是水均衡子上为多边形的立体柱体，就是水均衡子区。区。m+1mm-1(i,m)(j,m)(k,
30、m)3 3、多棱柱体三维流有限差分法、多棱柱体三维流有限差分法o 考虑以节点（考虑以节点（i i,m m）为中）为中心的多面柱体均衡子区，心的多面柱体均衡子区，下面根据达西定律和水下面根据达西定律和水均衡原理建立其差分方均衡原理建立其差分方程。如图所示，单位时程。如图所示，单位时间内通过单元中两断面间内通过单元中两断面pbpb、bq bq 断面流入到均衡断面流入到均衡区的水量为区的水量为:ijkpqacbxyo bqikhhTpbijhhTQnminmkmiknminmjmije1,1,1,1,(i,m)3 3、多棱柱体三维流有限差分法、多棱柱体三维流有限差分法,1,1,1,1,()44j m
31、j mi mi mij mij mzzzzTK,1,1,1,1,()44i mi mk mk mik mik mzzzzTKenminmkmiknminmjmijbqikhhTpbijhhTQ1,1,1,1,侧eQmikmijTT,式中：式中：所对应的柱面流入均衡区内的水量；所对应的柱面流入均衡区内的水量；为第为第m层流段层流段ij和和ik单元的导水系数，定义为单元的导水系数，定义为bqpb,从第从第e号三角形通过两线段号三角形通过两线段bqpbikij,流段的长度。流段的长度。类似对格点类似对格点 i 周围所有三角形作上述计算，并求和得到周围所有三角形作上述计算，并求和得到从侧面流入到均衡区
32、内的总水量从侧面流入到均衡区内的总水量侧QenminmkmiknminmjmijbqikhhTpbijhhTQ1,1,1,1,侧为为另外，容易求得通过均衡区顶底面流入到均衡区内的水量为另外，容易求得通过均衡区顶底面流入到均衡区内的水量为 1111,1,1,1/2,1/2,1,1,nnnni mi mi mi mViz miiz mii mi mi mi mhhhhQKAKAzzzz式中：式中：iA,1/2iz mK,1/2iz mKmiz,第第i个均衡区的平面投影区域面积；个均衡区的平面投影区域面积；第第i i个平面格点对应第个平面格点对应第m层和层和m-1-1层之间层之间z方向的平均渗透系数
33、；方向的平均渗透系数；第第i个平面格点对应第个平面格点对应第m层和层和m-1层之间层之间z方向的平均渗透系数；方向的平均渗透系数；第第i个平面格点对应第个平面格点对应第 m 层和层和 m-1层之间层之间z方向的平均渗透系数；方向的平均渗透系数；,111,1,1/2,1,i mnnnni mi mi mi miz miwi mii mi mhhhhKAQAzzt),.,2,1,.,2,1(iMmNi miwQ,若记若记 为第为第i个平面格点对应第个平面格点对应第m层的均衡子区的源汇项，层的均衡子区的源汇项，则由水均衡原理可获得三维地下水差分方程：则由水均衡原理可获得三维地下水差分方程：式中：式中
34、：N研究区平面网格结点总数，研究区平面网格结点总数，Mi第第i个平面结点对应的垂向分层数。个平面结点对应的垂向分层数。比储水系数比储水系数111111,1,1/2,1nnnnnnj mi mk mi mi mi mij mik miz miei mi mhhhhhhTpbTbqKAzzijik=Vi 对于非稳定流动问题，由初始时刻的水头分布开始，利用上时段的水头值递推计算以后各个时段的水头值，因此公式中为已知，为待求水头，将方程整理得：11,1/21,1/21,1,1,1,1,nni mini minij mj mik mk mi mi meei mi mi mi mKAKApbbqThThh
35、hzzzzijik1,1,2/1,1,2/1,nmiimimimiimimimiimiemikmijhtAzzAKzzAKikbqhTijpbT,i mi minwi mAQht),.,2,1,.,2,1(iMmNi-4、潜水（无压）水流有限差分法o 潜水含水层平面二维流()()xxyydHHHKHZKHZwxxyyt()()xxzzsHHHkkwxxzzto 潜水含水层剖面二维流o 对于剖面二维流，完全可以作为三维流的一种特例对于剖面二维流，完全可以作为三维流的一种特例处理，处理，因此这里只讨论平面二维流情形因此这里只讨论平面二维流情形。4.1显-隐式法 潜水含水层与承压含水层不同，随着潜水
36、面的升降，潜水含水层与承压含水层不同，随着潜水面的升降，会引起导水系数的变化，所以潜水含水层的导水系数是时会引起导水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓凝固系数法就是将某时段的潜水含水层的间的函数。所谓凝固系数法就是将某时段的潜水含水层的导水系数利用该时段初的水头值计算，即令导水系数利用该时段初的水头值计算，即令)(ZHKTxxxx)(ZHKTyyyy)()(yHTyxHTxtHyyxxd从而，潜水流方程可写成从而，潜水流方程可写成4、潜水（无压）水流有限差分法3.7.2 全隐式法)()()(ZHKTmxxmxx)()()(ZHKTmyymyyWyHTyxHTxtHmmyym
37、mxxmd)()()()()()()(潜水含水层与承压含水层不同，随着潜水面的升降，会引起导潜水含水层与承压含水层不同，随着潜水面的升降，会引起导水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓双重迭代法就是利用迭代法求解非线性潜水流微分方程。在第双重迭代法就是利用迭代法求解非线性潜水流微分方程。在第m m迭迭代步，潜水含水层的导水系数利用该迭代步的水头值计算，即令代步，潜水含水层的导水系数利用该迭代步的水头值计算，即令求解方程求解方程得到得到H H（m m），直到收敛为止。，直到收敛为止。4、潜水（无压）水流有限差分法3-8-1
38、 越流、入渗和抽水井等问题的处理越流、入渗和抽水井等问题的处理o如果考虑垂直渗流项（即源汇项），则二维承压流动微分方如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成程可写成o建立差分方程时，在结点处应加上这一项，它可具体表示为建立差分方程时，在结点处应加上这一项，它可具体表示为式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水层水头；为主含水层水头；为入渗强度或
39、蒸发强度；为抽水或注水流量。水或注水流量。5、特殊问题的处理：井孔水头校正2222HHHTTxy,()i juzi ji ji ji jzQKhhWMx y 井水位校正对比图井水位校正对比图 图图3-53-5（a a）大步长网格）大步长网格 图图3-53-5（b b）加密网格）加密网格有限差分法计算井水位的校正有限差分法计算井水位的校正 Q0 为注水为注水；Q0 为抽水为抽水 x=50mx=50mx=10m)2(ln2,wnjinwrxTQhh5、特殊问题的处理：、特殊问题的处理：不规则边界问题不规则边界问题o当研究区的几何形状属于简单形当研究区的几何形状属于简单形式（如矩形渗流区）时，差分网
40、式（如矩形渗流区）时，差分网格的划分往往将结点设在边界上。格的划分往往将结点设在边界上。o然而，对于实际问题来说，边界然而，对于实际问题来说，边界通常不是那么规则，边界的某些通常不是那么规则，边界的某些部分，甚至大部分不能与结点重部分，甚至大部分不能与结点重合，我们称这种边界是不规则的。合，我们称这种边界是不规则的。o关于不规则边界问题，直接取最关于不规则边界问题，直接取最靠近边界的曲折格线为近似边界靠近边界的曲折格线为近似边界越流、入渗和抽水井等问题的处理越流、入渗和抽水井等问题的处理l如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程如果考虑垂直渗流项（即源汇项），则二维承压流动微分方
41、程可写成可写成l建立差分方程时，在结点处应加上这一项，它可具体表示为建立差分方程时，在结点处应加上这一项，它可具体表示为式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水流量。水流量。源汇项的处理源汇项的处理:井孔水头校正井孔水头校正2222HHHTTxy,()i juzi ji ji ji jzQKhhWMx
42、 y 6、差分方程求解差分方程求解l一维显式差分格式 )(2)(2 )(2)(2111212212212331324322122321tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnnnnnnnnn网格个数为网格个数为ni直接求解直接求解6、差分方程求解差分方程求解 一维隐式差分格式一维隐式差分格式 )(2)(2 )(2)(2111211112212211121331321413122122131211tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnn
43、nnnnnnn网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW7、差分方程的收敛性和稳定性差分方程的收敛性和稳定性l截断误差截断误差：用差商代替微商时，地下水流动方程产生用差商代替微商时，地下水流动方程产生的误差为截断误差。的误差为截断误差。l收敛性收敛性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。精确解。则称该差分格式是收敛的。l稳定性稳定性：如果在求解差分方程过程中，某时间
44、步引入如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。大，则称该差分格式是稳定的。102一维一维显示格式显示格式的收敛条件和稳定条件是：的收敛条件和稳定条件是：2()()OtOx8、算例：、算例：显式有限差格式显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。各向同性。应用实例：河间地块承压水流模型应用实例：河间地块承压水流模型步骤：步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3
45、）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 8、算例：、算例：显式有限差格式（续显式有限差格式（续1）建立数学模型建立数学模型0,0tLxLx 022eHHKMtx)(00 xHHt)(),(210tHtHLXx（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取从而量化各变量的取值范
46、围。本例，取x-x-轴原点位于左端河，右侧轴原点位于左端河，右侧为正向，设两河流间距为为正向，设两河流间距为L.L.（3 3）数学模型）数学模型0t8、算例：、算例：显式有限差格式（续显式有限差格式（续2）将（0L）分成 N 等份，Nlx 1)1)网格剖分：网格剖分：xixittntn取时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、44）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格系统中任意取一点：在网格系统中任意取一点),(nitx),(txH),(nitx设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）22(
47、,)nneiniiHHKMx ttx8、算例：、算例：显式有限差格式（续显式有限差格式（续3）22(,)nneiniiHHKMx ttx1()nnniiiHHHOttt2211222()()nnnniiiiHHHHOxxx用差商代替微商用差商代替微商：将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记niHnih为为11122(,)()nnnnniiiiiinehhhhhKMx txt2)()(xOtO显然该式具有截断误差得到得到8、算例：、算例：显式有限差格式显式有限差格式（续4）111(12)(,)nnnniiiiinethhhhx t引入无量纲变量引入无量纲变量：将该式子代
48、入得到：将该式子代入得到：（i=1，2，3，.N-1），（n=1，2，3，.）2)(xtKMe11122(,)()nnnnniiiiiinehhhhhKMx txt8、算例：、算例：显式有限差格式显式有限差格式（续5）3）显示差分方程的求解）显示差分方程的求解l计算各结点初始时刻水头值计算各结点初始时刻水头值l利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值l利用边界条件计算边界结点水头值利用边界条件计算边界结点水头值l重复重复2 2、3 3步，直到计算出拟计算的各个时刻的水头步，直到计算出拟计算的各个时刻的水头值值)(00 xHHt8、算例：、算例：显式有限差格式（续
49、显式有限差格式（续6）22eHHKMtx)(00 xHHt)(),(210tHtHLXx004.0em10)(m,20)(m,20 m/d,821ttMK在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200200米，时间步长为米，时间步长为0.250.25天，分别计算各节点天，分别计算各节点各各时刻的时刻的水头值。水头值。0 x20m,0 x10m,)(0 xH8、算例：、算例：显式有限差格式显式有限差格式（续7）4/1)200(*004.025.0*20*8)(22xtKMeTime/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=100
50、0 m02010101010100.252012.5101010100.502013.7510.6251010100.752014.531 11.250 10.156 10101.002015.078 11.797 10.391 10.039 101.252015.488 12.266 10.654 10.117 10nieninininithhhh,111/)21()2(4141)4121(4111111nininininininihhhhhhh8、算例：、算例：显式有限差格式（续显式有限差格式（续8）Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=10
