第三讲-地下水数值模拟原理及建模方法和步骤课件.ppt

地下水数值模拟原理及地下水数值模拟原理及建模方法和步骤建模方法和步骤 内容提要内容提要绪论绪论一、地下水流有限差分法原理一、地下水流有限差分法原理二、地下水数值模型建模步骤二、地下水数值模型建模步骤三、建模所需要的基本资料三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流五、地下水系统热对流-扩散有限差分法原理扩散有限差分法原理绪论绪论地下水数值方法在水文地质学中的位置地下水数值方法在水文地质学中的位置o已学课程已学课程n水文地质学基础水文地质学基础n地下水动力学地下水动力学n水文地球化学（环境水文地质学）水文地球化学（环境水文地质学）n同位素同位素n地下水溶质运移：示踪剂或污染物在含水层中地下水溶质运移：示踪剂或污染物在含水层中的迁移机理及数学模型和求解方法的迁移机理及数学模型和求解方法o地下水数值模拟地下水数值模拟绪论绪论地下水数值方法在水文地质学中的位置地下水数值方法在水文地质学中的位置地下水动力学主要内容地下水动力学主要内容o连续性原理、达西定律、水均衡原理连续性原理、达西定律、水均衡原理o地下水流基本方程地下水流基本方程o几类特殊水文地质问题数学模型及解析解几类特殊水文地质问题数学模型及解析解n地下水向沟渠河中的流动地下水向沟渠河中的流动n园岛模型园岛模型n泰斯模型泰斯模型n有越流的不稳定井流（有越流的不稳定井流（Hantush and Jacob）n无越流的潜水含水层不稳定井流（无越流的潜水含水层不稳定井流（Neuman）o抽水试验及反求参数抽水试验及反求参数o数值方法可以作为地下水动力学课程内容之一、或补充、或延续数值方法可以作为地下水动力学课程内容之一、或补充、或延续绪论绪论数值解与解析解数值解与解析解o地下水动力学中所得到的解是解析解地下水动力学中所得到的解是解析解o解析解的特点解析解的特点n可以利用解析解公式计算出研究区任意空间点可以利用解析解公式计算出研究区任意空间点和时间点处的水头值和时间点处的水头值n对实际模型的概化对实际模型的概化n适用于一些特殊简单问题适用于一些特殊简单问题n公式复杂，需要借助计算机求解公式复杂，需要借助计算机求解o因此，其应用受到很大限制因此，其应用受到很大限制绪论绪论数值解与解析解数值解与解析解o数值方法是地下水动力学的完善和补充或延续数值方法是地下水动力学的完善和补充或延续o数值解的特点：数值解的特点：n只是求出研究区某些空间点和某些时间点处的水头值只是求出研究区某些空间点和某些时间点处的水头值n适用于所有的问题适用于所有的问题n具备水文地质基础和线性代数知识具备水文地质基础和线性代数知识n已有数值模拟专门软件（或自己编程）已有数值模拟专门软件（或自己编程）n需要有高性能计算机需要有高性能计算机n对实际问题的刻画比较精确对实际问题的刻画比较精确o因此，其应用非常广泛因此，其应用非常广泛数值方法很多，但是最简单实用的是有限差分法：数值方法很多，但是最简单实用的是有限差分法：绪论绪论数值方法数值方法 有限差分法有限差分法有限单元法有限单元法积分有限差分法积分有限差分法半解析半数值法半解析半数值法边界元法边界元法有限体积法有限体积法只讲有限差分法只讲有限差分法一、有限差分法的基本原理一、有限差分法的基本原理有两种方法建立差分方程有两种方法建立差分方程方法一方法一以地下水流基本微分方程及以地下水流基本微分方程及其定解条件为基础，其定解条件为基础，在渗流区剖分在渗流区剖分基础上，用差商代替微商，将地下水基础上，用差商代替微商，将地下水流微分方程的求解转化为差分方程流微分方程的求解转化为差分方程（代数方程）求解。（代数方程）求解。适用于二维矩形适用于二维矩形网格剖分、三维长方体网格剖分网格剖分、三维长方体网格剖分。方法二在渗流区剖分的基础上，直方法二在渗流区剖分的基础上，直接由达西定律和水均衡原理，建立各接由达西定律和水均衡原理，建立各个均衡区的水均衡方程，从而得到差个均衡区的水均衡方程，从而得到差分方程分方程。适用于矩形网格、三角形网适用于矩形网格、三角形网格。格。矩形网格矩形网格多边形网格多边形网格1、网格划分的基本类型、网格划分的基本类型o（1）先划格线，格点位）先划格线，格点位于网格中心于网格中心均均衡衡网网格格节节点点网网格格o（2）先规定格点位置，）先规定格点位置，再垂直平分两相邻结点的连再垂直平分两相邻结点的连线作格线，形成的网格即为线作格线，形成的网格即为水均衡区水均衡区方法一：方法一：差商代替微商差商代替微商MODFLOW网格系统网格系统方法一：方法一：差商代替微商差商代替微商导导数数的的有有限限差差商商近近似似导数的定义导数的定义 当当非常小的时候，有非常小的时候，有 上式右端项即为上式右端项即为f(x)(x)在在x0 0处的差商。处的差商。这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。方法一：方法一：差商代替微商差商代替微商(2)有限差分方程建立有限差分方程建立已知泰勒公式 由A得：AB 由B 得：称称 为为f(x)在在x0处的处的一阶前向差商，一阶前向差商，为为截断误差截断误差。称称 为为f(x)在在x0处的处的一阶后向差商，一阶后向差商，为为截断误差截断误差。方法一方法一 由A-B可以得：由A+B可以得：AB称称 为为f(x)在在x0处的处的一阶中心差商，一阶中心差商，为为截断误差截断误差。称称 为为f(x)在在x0处的处的二阶二阶中心中心差商，差商，为为截断误差截断误差。方法一方法一o对于偏导数（偏微商），类似可以得到相应的差商：对于偏导数（偏微商），类似可以得到相应的差商：方法一方法一(2)有限差分方程建立有限差分方程建立(续续)一维控制方程差分格式一维控制方程差分格式显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一控制方程控制方程网格剖分网格剖分nx个个二维控制方程差分格式二维控制方程差分格式显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一控制方程控制方程网格剖分网格剖分nx个个达西定律达西定律：水均衡原理水均衡原理：对某一研究对象，流入流入 流出流出体系内质量（或水量）变化量体系内质量（或水量）变化量研究对象可以是大区域的，也可以是微分单元体大区域的水均衡计算经常用于区域的水资源评价单位时间通过单位面积过水断面的水单位时间通过单位面积过水断面的水量与断面处的水力梯度成正比。量与断面处的水力梯度成正比。方法二：直接由达西定律和水均衡原理建立方法二：直接由达西定律和水均衡原理建立差分方程差分方程 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）+H 可取可取tn+1或或 tn时时刻的值刻的值节点（节点（i,j）的均衡区）的均衡区Aij:节点节点(i,j)均衡区的面积均衡区的面积(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例以二维承压水流为例 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）H取tn时刻的值时刻的值，有：+节点（节点（i,j）均）均衡区衡区Aij:节点节点(i,j)均衡区的面积均衡区的面积(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续1）(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy则得到显式格式：则得到显式格式：不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续2）+如果如果 H 取取tn+1时刻的值时刻的值节点（节点（i,j）的均衡区）的均衡区(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续3）记：记：得到隐式格式：得到隐式格式：(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1)xy不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续4）表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）当网格是等距时，即+H 可取可取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）的均衡区）的均衡区两边除以 ，得到 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）H 取取tn+1时刻的值，得时刻的值，得到隐式格式到隐式格式不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续6）表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）两边除以+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 表示：当水头下降一个单位时，从单位底面积饱水空隙介质中表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）释放的水量（体积）如果（i,j)是左边界节点，有+H 可取可取tn+1或或 tn时时刻的值刻的值节点（节点（i,j）的均衡区）的均衡区qij:节点节点(i,j)处的边界单宽流量。处的边界单宽流量。不等距矩形网格有限差分不等距矩形网格有限差分法：法：以二维承压水流为例（续以二维承压水流为例（续8）三维有限差分格式（三维有限差分格式（显示显示）三维有限差分格式（三维有限差分格式（隐式隐式）首先将渗流区划分成若干个辅助小三角形（图36中虚线所示）。划分时应注意：三角形的任一内角不得大于90，三条边的长度尽可能接近；三角形的顶点不能落在另外某个三角形的边上；相邻三角形变化不太大；应考虑水文地质条件，尽量使抽水井、观测孔（特别是用于拟合的观测孔）位于三角形的顶点处，河流、断层等位于三角形边上。1渗流区的剖分渗流区的剖分任意多边形网格系统任意多边形网格系统三角形的顶点称为格点，三角形的顶点称为格点，分为内格点和边界格点分为内格点和边界格点（包括第一类边界格点（包括第一类边界格点和第二类边界格点）。和第二类边界格点）。2、任意多边形网格有限差分法、任意多边形网格有限差分法一、有限差分法的基本原理一、有限差分法的基本原理2多边形均衡区网格的形成方法是：对某个格点i，分别作出格点i与其相邻格点的连线的垂直平分线，这些平分线所围成的多边形区域即为格点 i 的多边形均衡网格子区对每个格点都这样做，则形成多边形均衡网格系统（图36中实线所示）典型多边形子区典型多边形子区2、任意多边形网格有限差分法、任意多边形网格有限差分法 典型多边形子区典型多边形子区i=1,2,N隐式格式隐式格式水均衡方程为：水均衡方程为：2、任意多边形网格有限差分法、任意多边形网格有限差分法 典型多边形子区典型多边形子区i=1,2,N显式格式显式格式2、任意多边形网格有限差分法、任意多边形网格有限差分法o（1 1）区域剖分）区域剖分 首先，在平面上将研究区剖分成若干个辅助小三角形，首先，在平面上将研究区剖分成若干个辅助小三角形，然后以这些三角形为底面，将区域剖分成一些垂直的柱形体，然后以这些三角形为底面，将区域剖分成一些垂直的柱形体，再根据地层岩性及构造特征，用一些平面将柱形体剖分成若再根据地层岩性及构造特征，用一些平面将柱形体剖分成若干短柱体（图干短柱体（图3.13.1）。）。关于平面三角形的剖分依据下述原则：关于平面三角形的剖分依据下述原则：1 1）三角形的任）三角形的任一内角不得大于一内角不得大于9090，三条边的长度尽可能接近；，三条边的长度尽可能接近；2 2）三角形）三角形的顶点不能落在某个三角形的边上；的顶点不能落在某个三角形的边上；3 3）相邻三角形变化不太）相邻三角形变化不太大。大。3 3、多棱柱体三维流有限差分法、多棱柱体三维流有限差分法o（2 2）均衡区网格的形成）均衡区网格的形成 三角形的顶点称为格点，分为内格三角形的顶点称为格点，分为内格点和边界格点（包括第一类边界格点和点和边界格点（包括第一类边界格点和第二类边界格点）。第二类边界格点）。各结点的均衡区形成方法是：平面各结点的均衡区形成方法是：平面上，以格点为中心，与之相连的边的垂上，以格点为中心，与之相连的边的垂直平分线段组成一个多边形；垂向上，直平分线段组成一个多边形；垂向上，以该层与上下两层的两个中间层的平面以该层与上下两层的两个中间层的平面分别为上下底面，这样得到一个以平面分别为上下底面，这样得到一个以平面上为多边形的立体柱体，就是水均衡子上为多边形的立体柱体，就是水均衡子区。区。3 3、多棱柱体三维流有限差分法、多棱柱体三维流有限差分法o考虑以节点（考虑以节点（i i,m m）为中）为中心的多面柱体均衡子区，心的多面柱体均衡子区，下面根据达西定律和水下面根据达西定律和水均衡原理建立其差分方均衡原理建立其差分方程。如图所示，单位时程。如图所示，单位时间内通过单元中两断面间内通过单元中两断面pbpb、bq bq 断面流入到均衡断面流入到均衡区的水量为区的水量为:(i,m)3 3、多棱柱体三维流有限差分法、多棱柱体三维流有限差分法式中：式中：所对应的柱面流入均衡区内的水量；所对应的柱面流入均衡区内的水量；为第为第m层流段层流段ij和和ik单元的导水系数，定义为单元的导水系数，定义为从第从第e号三角形通过两线段号三角形通过两线段流段的长度。流段的长度。类似对格点类似对格点 i 周围所有三角形作上述计算，并求和得到周围所有三角形作上述计算，并求和得到从侧面流入到均衡区内的总水量从侧面流入到均衡区内的总水量为为另外，容易求得通过均衡区顶底面流入到均衡区内的水量为另外，容易求得通过均衡区顶底面流入到均衡区内的水量为 式中：式中：第第i个均衡区的平面投影区域面积；个均衡区的平面投影区域面积；第第i i个平面格点对应第个平面格点对应第m层和层和m-1-1层之间层之间z方向的平均渗透系数；方向的平均渗透系数；第第i个平面格点对应第个平面格点对应第m层和层和m-1层之间层之间z方向的平均渗透系数；方向的平均渗透系数；第第i个平面格点对应第个平面格点对应第 m 层和层和 m-1层之间层之间z方向的平均渗透系数；方向的平均渗透系数；若记若记 为第为第i个平面格点对应第个平面格点对应第m层的均衡子区的源汇项，层的均衡子区的源汇项，则由水均衡原理可获得三维地下水差分方程：则由水均衡原理可获得三维地下水差分方程：式中：式中：N研究区平面网格结点总数，研究区平面网格结点总数，Mi第第i个平面结点对应的垂向分层数。个平面结点对应的垂向分层数。比储水系数比储水系数=Vi对于非稳定流动问题，由初始时刻的水头分布开始，利用上时段的水头值递推计算以后各个时段的水头值，因此公式中为已知，为待求水头，将方程整理得：-4、潜水（无压）水流有限差分法o潜水含水层平面二维流o潜水含水层剖面二维流o对于剖面二维流，完全可以作为三维流的一种特例对于剖面二维流，完全可以作为三维流的一种特例处理，处理，因此这里只讨论平面二维流情形因此这里只讨论平面二维流情形。4.1显-隐式法 潜水含水层与承压含水层不同，随着潜水面的升降，潜水含水层与承压含水层不同，随着潜水面的升降，会引起导水系数的变化，所以潜水含水层的导水系数是时会引起导水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓凝固系数法就是将某时段的潜水含水层的间的函数。所谓凝固系数法就是将某时段的潜水含水层的导水系数利用该时段初的水头值计算，即令导水系数利用该时段初的水头值计算，即令从而，潜水流方程可写成从而，潜水流方程可写成4、潜水（无压）水流有限差分法3.7.2 全隐式法 潜水含水层与承压含水层不同，随着潜水面的升降，会引起导潜水含水层与承压含水层不同，随着潜水面的升降，会引起导水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓双水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓双重迭代法就是利用迭代法求解非线性潜水流微分方程。在第重迭代法就是利用迭代法求解非线性潜水流微分方程。在第m m迭代迭代步，潜水含水层的导水系数利用该迭代步的水头值计算，即令步，潜水含水层的导水系数利用该迭代步的水头值计算，即令求解方程求解方程得到得到H H（m m），直到收敛为止。，直到收敛为止。4、潜水（无压）水流有限差分法3-8-1 越流、入渗和抽水井等问题的处理越流、入渗和抽水井等问题的处理o如果考虑垂直渗流项（即源汇项），则二维承压流动微分方如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成程可写成o建立差分方程时，在结点处应加上这一项，它可具体表示为建立差分方程时，在结点处应加上这一项，它可具体表示为式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水流量。水或注水流量。5、特殊问题的处理：井孔水头校正井水位校正对比图井水位校正对比图 图图3-53-5（a a）大步）大步长长网格网格 图图3-53-5（b b）加密网格）加密网格有限差分法计算井水位的校正有限差分法计算井水位的校正 Q0 为注水为注水；Q0 为抽水为抽水 x=50mx=50mx=10m5、特殊问题的处理：、特殊问题的处理：不规则边界问题不规则边界问题o当研究区的几何形状属于简单形当研究区的几何形状属于简单形式（如矩形渗流区）时，差分网式（如矩形渗流区）时，差分网格的划分往往将结点设在边界上。格的划分往往将结点设在边界上。o然而，对于实际问题来说，边界然而，对于实际问题来说，边界通常不是那么规则，边界的某些通常不是那么规则，边界的某些部分，甚至大部分不能与结点重部分，甚至大部分不能与结点重合，我们称这种边界是不规则的。合，我们称这种边界是不规则的。o关于不规则边界问题，直接取最关于不规则边界问题，直接取最靠近边界的曲折格线为近似边界靠近边界的曲折格线为近似边界越流、入渗和抽水井等问题的处理越流、入渗和抽水井等问题的处理l如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成可写成l建立差分方程时，在结点处应加上这一项，它可具体表示为建立差分方程时，在结点处应加上这一项，它可具体表示为式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水流量。流量。源汇项的处理源汇项的处理:井孔水头校正井孔水头校正6、差分方程求解差分方程求解l一维显式差分格式网格个数为网格个数为ni直接求解直接求解6、差分方程求解差分方程求解一维隐式差分格式一维隐式差分格式网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW7、差分方程的收敛性和稳定性差分方程的收敛性和稳定性l截断误差截断误差：用差商代替微商时，地下水流动方程产生用差商代替微商时，地下水流动方程产生的误差为截断误差。的误差为截断误差。l收敛性收敛性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。精确解。则称该差分格式是收敛的。l稳定性稳定性：如果在求解差分方程过程中，某时间步引入如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。大，则称该差分格式是稳定的。一维一维显示格式显示格式的收敛条件和稳定条件是：的收敛条件和稳定条件是：8、算例：、算例：显式有限差格式显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。各向同性。应用实例：河间地块承压水流模型应用实例：河间地块承压水流模型步骤：步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 8、算例：、算例：显式有限差格式（续显式有限差格式（续1）建立数学模型建立数学模型（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取从而量化各变量的取值范围。本例，取x-x-轴原点位于左端河，右侧轴原点位于左端河，右侧为正向，设两河流间距为为正向，设两河流间距为L.L.（3 3）数学模型）数学模型8、算例：、算例：显式有限差格式（续显式有限差格式（续2）将（0L）分成 N 等份，1)1)网格剖分：网格剖分：取时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、4 4）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格系统中任意取一点：在网格系统中任意取一点设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）8、算例：、算例：显式有限差格式（续显式有限差格式（续3）用差商代替微商用差商代替微商：将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记为为显然该式具有截断误差得到得到8、算例：、算例：显式有限差格式显式有限差格式（续4）引入无量纲变量引入无量纲变量：将该式子代入得到：将该式子代入得到：（i=1，2，3，.N-1），（n=1，2，3，.）8、算例：、算例：显式有限差格式显式有限差格式（续5）3）显示差分方程的求解）显示差分方程的求解l计算各结点初始时刻水头值计算各结点初始时刻水头值l利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值l利用边界条件计算边界结点水头值利用边界条件计算边界结点水头值l重复重复2 2、3 3步，直到计算出拟计算的各个时刻的水头步，直到计算出拟计算的各个时刻的水头值值8、算例：、算例：显式有限差格式（续显式有限差格式（续6）在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200200米，时间步长为米，时间步长为0.250.25天，分别计算各节点天，分别计算各节点各各时刻的时刻的水头值。水头值。8、算例：、算例：显式有限差格式显式有限差格式（续7）Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 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 108、算例：、算例：显式有限差格式（续显式有限差格式（续8）Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252020101010100.502010201010100.7520300 2010101.0020-10 101.252010如果如果 t=1,t=1,则则8、算例：、算例：显式有限差格式（续显式有限差格式（续9）在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200米，时间步长为米，时间步长为0.25天，用隐式差分格式计算各节天，用隐式差分格式计算各节点个时刻的水头值。点个时刻的水头值。8、算例：、算例：隐式有限差格式隐式有限差格式在这个例子中，在这个例子中，解：隐式格式一般方程为解：隐式格式一般方程为于是有于是有根据初始条件得根据初始条件得根据边界条件得根据边界条件得由初始条件和边界条件由初始条件和边界条件由此解得由此解得t1t1时刻的水头值为时刻的水头值为在上述方程中取在上述方程中取 n=0,n=0,可以得到计算可以得到计算t1t1时刻水头值的方程时刻水头值的方程所以上述方程变成所以上述方程变成同理，可计算同理，可计算t2t2时刻的水头值时刻的水头值给定一组试探解给定一组试探解 ,分别代入上述三个分别代入上述三个方程，计算出方程，计算出l简单迭代法简单迭代法对于方程组对于方程组如果如果 是原方程组的解，则是原方程组的解，则 与与 非常接近。非常接近。8、算例：、算例：显式有限差格式（续显式有限差格式（续2）通常给定一个很小的正数通常给定一个很小的正数，若，若在计算出第在计算出第K步的迭代值后，按下式计算第步的迭代值后，按下式计算第K+1步的迭代值，步的迭代值，则称则称 为原方程组的近似解。为原方程组的近似解。一般情形，不可能一次就找到方程组的近似解，通常要计算很多一般情形，不可能一次就找到方程组的近似解，通常要计算很多步才可能得到满足精度要求的近似解。步才可能得到满足精度要求的近似解。若若 则取则取 为原方程为原方程组的近似解。组的近似解。8、算例：、算例：显式有限差格式（续显式有限差格式（续3）简单迭代法的收敛条件：如果方程组的系数矩阵是对角占优矩阵，简单迭代法的收敛条件：如果方程组的系数矩阵是对角占优矩阵，则相应的迭代格式是收敛的。由收敛迭代格式计算出的向量序列则相应的迭代格式是收敛的。由收敛迭代格式计算出的向量序列 一定收敛到某个确定的向量一定收敛到某个确定的向量 ，该向量就是原方程组的精确，该向量就是原方程组的精确解。解。8、算例：、算例：显式有限差格式（续显式有限差格式（续1）迭代法迭代法-改进迭代法改进迭代法在计算出第在计算出第K步的迭代值后，按下式计算第步的迭代值后，按下式计算第K+1步的迭代值，步的迭代值，一般情形，对于一个收敛的迭代格式来说，第一般情形，对于一个收敛的迭代格式来说，第k+1步的迭代值步的迭代值会比第会比第K步的迭代值更接近精确解，因此，根据这一认识，可以步的迭代值更接近精确解，因此，根据这一认识，可以对简单迭代格式进行改进。对简单迭代格式进行改进。若若 则取则取 为原方为原方程组的近似解。程组的近似解。迭代法迭代法-改进迭代法（续）改进迭代法（续）通常给定一个很小的正数，若通常给定一个很小的正数，若在计算出第在计算出第K步的迭代值后，按下式计算第步的迭代值后，按下式计算第K+1步的迭代值，步的迭代值，则称则称 为原方程组的近似解。为原方程组的近似解。一般情形，不可能一次就找到方程组的近似解，通常要计算很多一般情形，不可能一次就找到方程组的近似解，通常要计算很多步才可能得到满足精度要求的近似解。步才可能得到满足精度要求的近似解。若若 则取则取 为原为原方程组的近似解。方程组的近似解。内容提要内容提要绪论绪论一、地下水流有限差分法原理一、地下水流有限差分法原理二、地下水数值模型建模步骤二、地下水数值模型建模步骤三、建模所需要的基本资料三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流五、地下水系统热对流-扩散有限差分法原理扩散有限差分法原理二、地下水数值模型建模步骤二、地下水数值模型建模步骤模拟步骤模拟步骤o建立概念模型建立概念模型 o建立数学模型建立数学模型 o数值方法及软件数值方法及软件（编程）（编程）o参数准备以及计算参数准备以及计算 o模型校正与检验模型校正与检验 o参数敏感性分析参数敏感性分析 o预测模拟预测模拟软件软件一、概念模型（模型概化）一、概念模型（模型概化）根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿物、水文、气象、工农业利用情况等岩石矿物、水文、气象、工农业利用情况等模拟的区域模拟的区域:含水层类型含水层类型:潜水（无压）、承压、混合、多层潜水（无压）、承压、混合、多层维数维数:一维、二维、三维一维、二维、三维水流状态：稳定流水流状态：稳定流/非稳定流、饱和流非稳定流、饱和流/非饱和流非饱和流介质状况介质状况:均质和非均质均质和非均质/各向同性和各向异性各向同性和各向异性孔隙孔隙/裂隙裂隙/双重介质双重介质流体的密度差流体的密度差边界条件和初始条件边界条件和初始条件 必要时需进行一系列的室内试验与野外试验必要时需进行一系列的室内试验与野外试验,以获取有关参数以获取有关参数,如渗透系如渗透系数、弥散系数、分配系数、反应速率常数等。数、弥散系数、分配系数、反应速率常数等。收集研究区已有收集研究区已有水文地质资料水文地质资料概念模型概念模型边界概化边界概化内部结构概化内部结构概化完成模型概化图完成模型概化图地表水体地表水体水头边界水头边界流量边界流量边界断层接触边界断层接触边界隔水边界隔水边界水头边界水头边界流量边界流量边界岩体岩层接触边界岩体岩层接触边界天然分水岭天然分水岭据含水层组类型、据含水层组类型、结构、岩性结构、岩性含水介质含水介质地下水运动状态地下水运动状态水文地质参数水文地质参数确定层组的均质、非均质，确定层组的均质、非均质，各向同行、异性，稳定、非各向同行、异性，稳定、非稳定，潜水、承压水稳定，潜水、承压水空间分布空间分布导水、储水系数、主渗透方向导水、储水系数、主渗透方向目的层与相邻层关系目的层与相邻层关系均质、非均质，各向同、异性均质、非均质，各向同、异性层、紊流，二、三维层、紊流，二、三维时间概化时间概化渗透、储水系数，给水度及单位涌渗透、储水系数，给水度及单位涌水量，含水层分布规律，地下水流水量，含水层分布规律，地下水流场，水化学场、温度场的空间概化场，水化学场、温度场的空间概化平面图平面图剖面图剖面图研究范围、主要居民点、标志性地形、地貌、地表河流、研究范围、主要居民点、标志性地形、地貌、地表河流、湖泊、开采井、地下水天然露头、含水层的各类钻孔湖泊、开采井、地下水天然露头、含水层的各类钻孔地表地理要素、含水层结构、地质结构、地下水水位、地表地理要素、含水层结构、地质结构、地下水水位、各类源汇项及其性质各类源汇项及其性质（二）数学模型（二）数学模型 三维地下水流动问题控制方程三维地下水流动问题控制方程第二类边界条件第二类边界条件第一类边界条件第一类边界条件初始条件初始条件 绝大部分数学模型是无法用解析法求解的，数绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。值化就是将数学模型转化为可解的数值模型。三、数值方法及软件（或编程）三、数值方法及软件（或编程）有限差分法有限差分法有限单元法有限单元法积分有限差分法积分有限差分法半解析半数值法半解析半数值法边界元法边界元法有限差分法有限差分法:MOFLOW系列系列GMS 中中MODFLOWVisual MODFLOWProcessing MODFLOW有限单元法有限单元法:FEFLOW积分有限差分法积分有限差分法:TOUGH2，TOUGH REACT三、数值方法及软件（或编程）三、数值方法及软件（或编程）四、模型参数准备以及计算四、模型参数准备以及计算o含水层参数：含水层参数：渗透系数，渗透系数，弹性释水系数，孔隙度等弹性释水系数，孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、大气降水入渗系数（分区、数值）数值）n蒸发排泄系数蒸发排泄系数n地表水体水位、底面高程、地表水体水位、底面高程、底面岩性特征底面岩性特征n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状）人工开采（点状、面状）o边界条件边界条件o初始条件初始条件计算计算资料资料结果结果输入模块输入模块运行模块运行模块输出模块输出模块Visual MODFLOW参数不确定性参数不确定性钻孔太少，钻孔太少，地层资料少，地层资料少，钻孔多，含钻孔多，含水层结构会水层结构会发生变化发生变化四、模型校正与检验四、模型校正与检验（1）钻孔资料）钻孔资料（2）抽水试验）抽水试验参数不确定性参数不确定性泰斯井流抽水试验假设：泰斯井流抽水试验假设：（1）（2）（3）（4）（5）（6）（四）模型校正与检验（四）模型校正与检验（1）钻孔资料）钻孔资料（2）抽水试验）抽水试验（五）模拟：模型校正（参数识别（五）模拟：模型校正（参数识别）将模拟结果与实测结果比较将模拟结果与实测结果比较,进行参数调整进行参数调整,使模使模拟结果在给定的误差范围内与实测结果吻合。拟结果在给定的误差范围内与实测结果吻合。调参过程是一个复杂而辛苦的工作调参过程是一个复杂而辛苦的工作,所调整的参数所调整的参数必须符合模拟区的具体情况。必须符合模拟区的具体情况。人机交互与自动调参相结合。尽管自动调参程序人机交互与自动调参相结合。尽管自动调参程序(如如PEST),PEST),也不能代替人的工作。也不能代替人的工作。（五）模拟：模型检验（五）模拟：模型检验 模型验证是在模型校正的基础上模型验证是在模型校正的基础上,进一步进一步调整参数调整参数,使模拟结果与第二次实测结果吻合使模拟结果与第二次实测结果吻合,以进一步提高模型的置信度。以进一步提高模型的置信度。（六）灵敏度分析（六）灵敏度分析 校正后的模型受参数值的时空分布、边界校正后的模型受参数值的时空分布、边界条件、水流状态等条件、水流状态等不确定度不确定度的影响。灵敏度的影响。灵敏度分析就是为了确定不确定度对校正模型的影分析就是为了确定不确定度对校正模型的影响程度。响程度。（七）预测（七）预测 用校正的参数值进行预测用校正的参数值进行预测,预测时需估算未预测时需估算未来的水流状态。来的水流状态。后续检查与模型的重建（完善）后续检查与模型的重建（完善）后续检查在模拟研究结束数年后进行。收集新的野后续检查在模拟研究结束数年后进行。收集新的野外数据以确定预测结果是否正确。如果模拟结果精确外数据以确定预测结果是否正确。如果模拟结果精确,则该模型对该模拟区来说是有效的。由于场址的唯一则该模型对该模拟区来说是有效的。由于场址的唯一性性,故模型只对该模拟区有效。故模型只对该模拟区有效。后续检查应后续检查应在预测结束足够长的时间后进行在预测结束足够长的时间后进行,以以便有足够的时间发生明显的变化。便有足够的时间发生明显的变化。模型的再设计模型的再设计 一般来说一般