地下水运动的数学模型

1、第四章 地下水运动的数值模型解析解虽然具有精确可靠的特点，但采用解析解反映自然状态和复杂人类活动干扰下的地下水运动是相当困难的。因此，当含水层的条件严重偏离现有解析模型的简化假设时，人们通过数值模型来获得近似的地下水流场及演变趋势。第一节 地下水流数值方法概述地下水流的数学模型采用偏微分方程描述地下水流的时间和空间连续状态，而数值模型则是采用离散（非连续）时空模型中水头的分布与演变对数学模型进行近似描述。从精确数学模型到近似数值模型的转化，虽然会损失一些精度，但使复杂地下水流问题的分析得以通过机械计算实现，而且误差也是可控的。把偏微分方程求解的数值方法引入到地下水流问题的求解始于20世纪70年

2、代，主要方法包括有限差分法、有限元法和边界元法，此后又发展了有限分析法、多重网格法和无网格法等。这些方法的共同特点是将模型空间及边界离散为由一系列的节点以及联系这些节点的单元（无网格法除外），含水层的水头在这些节点上定义，从而实现了水头分布空间连续函数向离散变量的转化，表示为 （4.1.1）式中；H为含水层的水头；x、y、z为空间坐标；p为数值模型的节点；M为节点的数目。与此同时，时间也被离散化为一系列具有先后顺序的时刻，不同时刻节点的水头分布方式也不同，每个时刻节点的水头分布结果总是由前一个时刻演变过来，即 （4.1.2）式中：k为时刻；t为时间；为时间步长。数值模型的核心部分，是在地下水流

3、数学模型的基础上确定节点在流场中的相互关系，并形成如下形式的方程组（4.1.3）式中：系数反映了节点P与节点L在流场中的关系，可称为关联系数；通常是由前一个时刻决定的常量。式（4.1.3）除包括地下水流控制方程的离散近似以外，还包括了边界条件的离散化处理。当k=0时，中包含了初始条件。如果模拟稳定流，则与时间无关，上标k可以略去。使用矩阵符号，式（4.1.3）可缩写为 （4.1.4）实际在处理模型时，一个节点通常只与相隔最近的若干节点有直接的关系，在式（4.1.3）中大多数关联系数为零，而且一般有对称性，即。因此，矩阵通常是对称的稀疏矩阵，降低了方程组求解的困难。尽管如此，区域地下水流数值模型

4、的节点数往往很大，而且在某些非线性条件下，本身是待求水头的函数，方程组的求解需要使用迭代法而耗费相当多的时间。不同类型的地下水流数值模型，最核心的差异在于形成式（4.1.3）的方式，关联系数具有不同的计算公式。原则上，对于完全相同的地下水流数学模型，以及完全相同的离散节点序列和离散时间序列，不同的数值模型应给出近乎同样的结果，只是计算精度略有差异。已有研究表明，某种类型的有限差分法与某种类型的有限元法存在一定的等价性，但更多的对比分析尚待进一步探讨。实际上，研究者在检验数值模型的可靠性时，总是与解析解的模型对比。当不同方法的数值模型用相同的解析解作验证时，他们之间的等价性就可以得到一定程度的证

5、明。然而，人们对数值模型的偏好并不取决于数值模型之间的等价性，而往往取决于某种数值方法的熟悉程度，或者基于某种数值方法专业软件的方便快捷。专业软件在地下水流数值模型的推广应用中发挥了重要作用，如有限差分模拟程序MODFLOW（McDonald and Harbaugh，1988）和HST3D（Kipp，1987）、有限元模拟程序FEMWATER（Lin et al，1997）和FEFLOW（Diersch and Kolditz，1998）等，国内有陈崇希等开发的有限差分模拟程序PGMS（陈崇希等，2007）。其中，以美国地质调查局发布的有限差分模拟程序MODFLOW最为著名，本章最后将介绍其

6、特点。第二节 地下水流有限差分模型有限差分法是地下水流数值模拟最早使用的一种方法。这种方法古典而简洁、物理意义明了、编程容易，因而被广泛流传。地下水流的有限差分模型已经从最初的一维均质模型发展为目前的三维非均质模型，并有专业软件加以实现。一、有限差分法的基本原理有限差分法是求解偏微分方程边值问题和初值问题的一种数值方法，其实质是把连续的模型空间离散化为规则或不规则的网格点，利用导数的差分近似形式代替偏微分方程形成差分方程组，通过求解方程组得到离散点的待求变量作为连续场的一种近似结果。模型空间的离散化从形态上主要分为规则网格与不规则网格。在二维平面空间，规则网格由一系列与坐标轴平行的直线组成，直

7、线交叉形成的单元（cell）为矩形；在三维空间，规则网格由一系列与坐标轴正交的平面组成，平面相互切割形成的块体（block）为长方体。不规则网格在平面上由一系列互不重叠的多边形填充模型空间而形成，一般通过分层建立三维空间网格。一维问题空间离散化很简单，就是在坐标轴上生成一系列具有一定间隔的点，或称节点（node），相邻两个节点之间为线段，待求变量在这些节点上的分布代表了在坐标轴上连续分布的近似状态。有限差分法中的“有限”是指网格中的节点、单元或块体数目是有限的，每个线段、单元或块体的尺度也是有限的。首先用一维问题来阐述导数的差分格式。设待求变量沿x轴的分布为二阶连续函数，现将x轴离散化为节点，

8、它们的待求变量为，则在节点处的一阶导数表示为 （4.2.1）式中：，为节点的间距。去掉式（4.2.1）中的Taylor展开式余项，就得到一阶导数前向差分格式。一阶导数还可以表示成后向差分格式为 （4.2.2）式中：。显然，第一个节点处的一阶导数只能表示为前向差分格式，而最后一个节点处的一阶导数只能表示为后向差分格式。有了一阶导数的差分格式之后，再次使用中心差分格式，就可以建立二阶导数的差分格式为 （4.2.3）其中：一阶导数在节点采用后向差分格式，在节点采用前向差分格式，则可以得到二阶导数的差分格式为 （4.2.4）如果相邻节点的间距都相等，即，则上式可简化为 （4.2.5）建立了一阶和二阶导

9、数的差分格式，就可以用它们来改写所研究问题的偏微分方程。下面用一个简单的例子来说明有限差分模型的建立方法。设有一边值问题的数学模型为 （4.2.6） （4.2.7） （4.2.8）在建立有限差分模型时，把模型空间0，2L离散化为三个节点。对第二个节点建立控制方程（4.2.6）的差分格式为 （4.2.9）第一个节点为边界节点，直接利用边界条件式（4.2.7）得 （4.2.10）第三个节点也是边界节点，把边界条件式（4.2.8）表示成差分格式为 （4.2.11）这样，由式（4.2.6）至式（4.2.8）描述的数学模型就转化为了由方程组式（4.2.9）至式（4.2.11）构成的差分模型。求解这个差分

10、方程组，得 （4.2.12）待求变量的差分模型结果就反映了数学模型精确解得特征。如果我们遇到的是一个初边值问题，例如，把上述关于的边值问题改写为如下关于的数学模型 （4.2.13） （4.2.14） （4.2.15） （4.2.16）其中，式（4.2.16）是初始条件。对于这种初边值问题，我们需要先解决关于时间的偏导数。把时间坐标轴离散化为时步，则第时刻的控制方程可以写为 （4.2.17）式中：为时间步长。时间偏导数取后向差分格式，方程左侧偏导数的自变量全部取第k时刻的数值，这种时间上的处理称为隐式差分格式。控制方程也可以处理为显式差分格式，即 （4.2.18）这样可以直接根据前一个时刻的变量

11、分布计算当前时刻的，为 （4.2.19）关于时间的偏导数问题解决以后，就可以把上述初始问题的求解转化为一系列离散时刻关于的边值问题的求解。例如，采用隐式差分格式，则的数学模型为 （4.2.20） （4.2.21） （4.2.22）利用边值问题的有限差分法，可以把上述数学模型转化为有限差分模型，得到各个节点上的待求变量，只要把前一时刻求解出来的结果fk-1代入到式（4.2.20）即可。数值模型求解的流程（图4.1）。k=0差分模式：k = M结束否是图4.1 初边值问题有限差分模型的求解流程示意图有限差分模型必须满足收敛性和稳定性才是有价值的。收敛性的含义是当离散网格的x、y、z、t趋于零时，有

12、限差分模型的离散解应趋于数学模型的精确解，即随着网格剖分尺度和时间步长的缩小，模型的误差也越来越小。稳定性则与初边值问题在时间上的推进有关，初边值问题在有限差分模型中是转化为各个时步的边值问题求解的，如果前一个时步差分方程求解带来的误差不会在以后的时步中积累，则差分模型是稳定的。否则模型不稳定，随着时步的推进，模型误差将不断积累增大，越来越偏离精确解。对地下水流的有限差分模型，人们已经证明（陈崇希等，1989）隐式差分格式是无条件收敛和稳定的，而显式差分格式虽然满足收敛性，但一般只有在时间步长较小、网格剖分尺度相对较大的时候才满足稳定性。因此，绝大多数地下水流有限差分模型采用隐式差分格式进行求

13、解。二、基于矩形网格的地下水流差分方程用矩形网格对模型空间进行离散，是地下水流数值模拟的常用方法，因为规则网格便于图形显示和数据处理。在矩形网格上建立二维地下水流的有限差分模型，既可以直接套用偏导数的差分格式，也可以根据水均衡原理进行推导。考虑给一个承压含水层建立矩形网格的差分模型，在矩形网格的每个单元形心处设置一个节点，如图4.2中的圆点，节点的离散网格编号为（i，j），其坐标为（xi，yj）。对于模型边界以内的节点，一般都有8个相邻节点，如图42中节点O（i，j）分别与节点18相邻。如果含水层的渗透性主轴方向与网格坐标轴的方向一致，则仅用节点l4就可以建立节点O所在矩形单元的地下水流差分方

14、程。图4.2 矩形网格差分模型的局部特征设含水层均质各向异性，导水系数为Tx和Ty，储水系数为S。在t=tk时刻，地下水的水头在节点上的分布为Hk，根据达西定律，从周围四个节点（节点l4）流向节点O的流量用差分格式可分别表示为而第k个时步内，节点O所在矩形单元内的地下水储存量增加值表示为 （4.2.23）这个方程中隐含的假定是：单元内的平均水头变化与节点0处的水头变化相等。根据水均衡原理，单元内水体积的增加量应与侧向流量在时间步长内的积累和其他外部补给量相互平衡，即 （4.2.24）式中：为其他外部补给流量。上式的展开式为（4.2.25）这就是对应节点0（i，j）的地下水流差分方程，它可以进一

15、步变形为 （4.2.26）其中：对于任何一个内部节点都可以写出形如式（4.2.25）的差分方程，然后运用边界条件写出边界节点的差分方程。所有节点的差分方程联立求解，即得到tk时刻水头分布的离散形式。如果含水层是非均质的，那么可以对每个节点所在的矩形单元设置不同的参数（导水系数T与储水系数S），例如，用分别表示节点O及相邻节点所在矩形单元的x方向和y方向的导水系数。那么节点之间的单位界面流量，可以用一个平均意义上的导水系数进行计算，以节点1流向节点O的流量为例，即 （4.2.27）其中：平均导水系数通常采用调和平均值 （4.2.28）依此类推其他节点与节点O之间的平均导水系数。考虑这种非均质参数

16、分布之后，式（4.2.26）中的各项系数更新为注意：其中的外部补给流量可能既与节点有关，又随时间变化。如果研究对象是潜水含水层，把上述差分方程中的储水系数S改为给水度Sy，至于导水系数T，则可以根据饱和厚度水头计算出每个节点所在单元的等效导水系数，以节点0为例，有 （4.2.29）式中：为含水层在节点0处的渗透系数；为含水层底板在节点0处的高程。这样就导致了一个新问题，在潜水含水层的差分方程（4.2.26）的中，系数C0C4等与待求的水头序列Hk有关，从而使式（4.2.26）变成非线性方程。解决这个问题的办法有以下几种。1）显式法：用前一个时刻的已知水头序列Hk-1计算潜水含水层的导水系数，在

17、时步内不再更新。当时间步长较大，水头变化明显时，这种方法精度较低。2）预测-校正法：首先用显式法求出tk-1/2时刻的水头序列Hk-1/2，时间步长设置为tk/2，然后计算出中间时刻的导水系数，重新计算系数C0C4，返回去由前一时刻的已知水头序列Hk-1求出待求水头序列号Hk，时间步长恢复为tk。3）迭代法：首先用显式法求出水头序列Hk，并计算出新的导水系数，重新计算系数C0C4，返回去再次求解水头序列Hk，如此反复，直到前后两次计算得到的水头几乎相同（由允许误差控制）为止。这种方法较为精确，但计算量大。基于矩形网格的地下水流三维模型，一般是进行分层处理，每层都是同样的矩形网格。网格单元为长方

18、体，每个长方体内设置节点，编号的形式为（i，j，l），其中的l表示层序号。在建立地下水流差分方程时，除了同层的4个相邻长方体单元之外，还需要考虑上下相邻的2个长方体单元，因此每个内部节点一般包含6个相邻节点。关下分层矩形网格的三维地下水流差分模型，在有关MODFLOW模型的章节中将具体阐述。三、基于多边形网格的地下水流差分方程基于矩形网格的有限差分法虽然简单，但是在处理不规则的模型边界和局部加密方面存在缺陷，会导致网格浪费。基于多边形网格的有限差分可以适应不规则的边界，并且能够在不显著增加整体节点数的情况下加密局部网格。多边形网格的形成，需要先建立模型空间的三角形剖分离散网格，作为辅助网格，再

19、用三角形单元各边的中垂线连接为各个节点所对应的多边形单元（图4.3）。图4.3 一个平面多边形网格的局部图（灰色线为辅助三角形网格）辅助三角形网格的剖分应注意（陈崇希等，1989）：三角形的任一内角应小于90，尽可能接近正三角形； 三角形的顶点不能落在另外某个三角形的边上，即节点之间的连线最多与2个三角形相连；节点的分布需要考虑地下水流动的控制条件，如河流、井点等。辅助三角形生成之后，再组建多边形网格，其方法类似于气象水文学中降雨量插值的泰森多边形法（Bear，1979），对每条三角形的边作垂直平分线，直到与其他的垂直平分线相交，这些垂直平分线连接成一个一个的多边形，每个多边形包围一个节点，作

20、为该节点的控制区。对具有多边形特征的节点控制区建立地下水流偏微分方程的差分格式，可以利用封闭域内体积分与边界积分之间的关系（高斯定律），其方式与有限体积法一致。因此，在某种程度上说，基于多边形网格的有限差分法属于有限体积法。不过，节点控制区的差分方程也可以在水均衡原理的基础上建立。设模型研究的对象是一个均质各向同性的承压含水层，从多边形网格的内部取一个节点控制区进行分析（图4.4）。图4.4 平面多边形网格的一个节点控制区节点0与N个节点相邻，对于第i个节点，节点0与节点i之间连线的长度为Li，而对应的垂直平分线的宽度为wi，则在节点i与节点0之间流向控制区的流量Qi可根据达西定律计算为 （4

21、.2.30）式中：T为含水层的导水系数；Hi和H0分别为节点i与节点0的水头。所有相邻节点与控制区之间都可以建立这样的关系，因此从控制区周围流向控制区内的地下水的总流量应为 （4.2.31）式（4.2.31）中把时间节点k也放进去了。在时间步长内，控制区内地下水的总体积的增加量为 （4.2.32）式中：S为含水层储水系数；A0为节点0的控制区面积，式（4.2.32）中隐含的假设是节点0处的水头变化与控制区内的平均水头变化相等。根据水均衡原理，控制区水体积的增加量应该与流入控制区的总水量相互平衡，即 （4.2.33）式中：为其他的外部补给流量。把式（4.2.33）展开得 （4.2.34）此式可以

22、进一步改写为 （4.2.35）这就是节点控制区的地下水流差分方程。对于多边形网格中的每个内部节点，都可以建立相同形式的差分方程，与边界节点的差分方程一起形成求解模型的方程组。每个相邻节点对应的wi/li数值，以及节点控制区的面积A0，可以根据辅助三角形网格的几何信息进行计算，具体方法见有关文献。利用多边形网格建立三维有限差分模型的方法，一般是建立分层的多边形网格以模拟含水层的构造特征图4.5（a）。每个节点除了和同一层的若干节点相邻之外，还与上层或下层相同平面位置的节点相邻，以图4.5（b）中的节点0为例，其同层相邻节点为节点1N，而节点u和节点l分别为上层相邻节点和下层相邻节点，它们与节点0

23、连线的中垂面和平面控制区共同形成节点0的三维控制体，这个控制体的水量均衡方程为 （4.2.36）式中：Qu和Qi分别为控制体顶面和底面流向控制体内的流量；Ss为含水层的储水率；V0为节点O控制体的体积。控制体顶底面的垂向流量根据达西定律计算为 （4.2.37）式中：为含水层的垂向渗透系数；，为节点u与节点0的高差；，为节点l与节点0之间的高差图4.5（b）。 图4.5 三维分层多边形网格及其节点控制区在式（4.2.36）中，计算侧向流量Qi时应使用水平渗透系数KH=Kx=Ky（仅考虑横观各向同性条件），每个节点控制区都可以计算出等效的导水系数，形如 （4.2.38）在计算侧向流量时使用调和平均

24、的导水系数，即 （4.2.39）把式（4.2.39）和式（4.2.37）代入式（4.2.36）可得三维网格地下水流的差分方程。四、地下水流控制条件的差分方程在建立了内部节点的差分方程之后，还必须建立地下水流控制条件的差分方程，包括边界条件和某些特殊的水文地质条件，其方法是建立具有不同边界类型的边界节点的差分方程，同时对某些具有特殊水文地质条件的内部节点改写差分方程。常见的控制条件包括以下类到。（1）已知水头边界如果节点在已知水头边界上，则直接把计算时刻的已知水头赋予该节点。（2）已知流量边界如果一个节点在已知流量或流速的边界上，则应把边界流量纳入差分方程（4.2.25）的外部补给项之中。设边界

25、上的单宽流量为q，沿着x坐标轴方向，则对于矩形网格图4.6（a）中的节点0，外部补给项计算为 （4.2.40）而差分方程（4.2.26）变为 （4.2.41） 式中：F已经包含了外部补给项.对于多边形网格图4.6（b）中的节点O，流向控制区的边界流量需要根据控制区外侧边在y轴上的投影长度L8计算，即 （4.2.42）节点1和节点4都属于边界节点。节点0的控制区对应的水流差分方程仍然保持式（4.2.35）。图4.6 已知流量边界上的节点（3） 越流控制条件如果要模拟第一类越流系统的地下水流，且只把承压含水层作为模拟层，则越流需要考虑在外部补给项中，即 （4.2.43）式中：和分别为弱透水层的垂向

26、渗透系数和厚度；是作为补给源的含水层的固定水头；为节点0所在矩形单元或控制区的面积。对于矩形多边形网格，越流量加入方程（4.2.25）中，但由于其中含有待求变量节点0的水头，差分方程（4.2.26）中的系数C0和常数项F需要修改为 （4.2.44） （4.2.45）式中：是除了越流量以外其他的外部补给流量。（4） 潜水面在三维有限差分模型中，具有潜水面的模拟层一般被作为顶层处理，其特殊之处在于导水系数T用地下水饱和厚度计算，而地下水储存量的变化部分用给水度SY计算，导致模拟层的节点差分方程为非线性方程，可用显示法、预测-校正法或迭代法求解。在实际模拟过程中，往往还有可能出现顶部模拟层被疏干、而

27、潜水面移动到下部模拟层的现象，需要进行技术上的处理。目前，不同的专业模型软件对此具有不同的处理方法，如MODFELOW模型中使用Wetting和Drying模块技术。（5） 井孔井孔在模拟含水层-井孔系统地下水流时，是最重要的控制条件之一。然而，常规的有限差分模型在处理抽水井和定降深井方面存在缺陷，需要进行校正，具体方法见下一小节。（6） 其他特殊条件地下水与河流、湖泊和泉点等的相互作用，降水入渗补给以及潜水的蒸发排泄等，都是需要在有限差分模型中处理的特殊条件。各种处理方法属于地下水流数值模拟的专门技术，反映了模拟者对这些特殊条件下地下水流物理机制的认识。五、地下水流差分模型的井孔校正常规的地

28、下水流数值模型一般仅仅把井孔的抽水流量以外部补给项Qa（抽水为负）的形式加入差分方程，虽然可以反映抽水井对周围节点的影响，但抽水井本身所在位置的水头结果不正确。另外，对于定降深井，如自流井，常规的做法是直接把井孔的已知水头设置到所在的节点上，但计算得到的井流量不正确。这是一种“算不准”现象，即已知流量的井孔算不准井孔水头，而已知水头的井孔算不准井孔流量。为了解决“算不准”问题，需要对常规的有限差分模型进行校正。关于抽水井的问题，Prickett（1967）、Peaceman（1978）和陈崇希等（1989）已经指出有限差分模型中井点单元（well-block）的计算水头不等于井孔中的水头，通过

29、引入井孔附近稳定流或拟稳定流的解析解，提出了校正公式。以正方形网格为例（图4.7），校正公式（陈崇希等，1989）为 （4.2.46）可以通过引入一个等效半径来改写这个公式，即 （4.2.47）其含义是：差分模型求解得到的井点单元水头，相当于距离井心为处的水头，并非实际的井孔水头。等效半径只与差分网格的尺寸有关，即 （4.2.48）由于井孔半径一般远小于模型的网格尺寸，模型计算水头往往远远低估了抽水井本身的降深，在绘制井孔周围等值线图时显示出过小的水力梯度。有了校正式（4.2.47）之后，就可以先用常规优先差分模型计算井点单元的水头，然后直接校正得到井孔的水头。在校正时，是必须提供的参数。图4

30、.7 正方形差分网格及抽水井附近水头分布图然而，上述校正公式是在正方形网格、无面源补给的情况下调用稳定井流解析解推导出来的，是否能够适应矩形网格、不规则网格以及有面源补给（越流等）的情况呢？定降深井孔是否也能直接用式（4.2.47）进行校正呢？在此，我们对承压含水层建立一种更加普遍的校正方法来解决这个问题。为了得到具有实用性的解析解，对井孔附近采取以下的简化假设：1) 井孔附近的地下水在水平方向上是径向流动；2) 井孔附近的含水层为均质各向同性，导水系数为T、储水系数为S；3) 井孔的滤管穿透了所研究的含水层（模拟层）。在上述假设的基础上，井孔附近的地下水流可以用下述方程为 （4.2.49）式

31、中：是一个由外部因素决定的面源补给项。引入一个因子改写上述方程，即 （4.2.50） （4.2.51）这个因子综合考虑了面源补给与含水层储存量的变化，并且我们再引入一个简化假设：因子在井孔附近与径向距离r无关。这个假设并不符合实际，目的是在当前条件下获得一个实用的近似解。给定流量为，水头为的井孔边界条件 （4.2.52）可以求出方程（4.2.50）满足因子假设的解析解为 （4.2.53）由于井孔半径一般远小于研究尺度下的径向距离r，上式可简化为 （4.2.54）对于矩形网格（图4.8），常规差分模型采用下式建立井点0所在矩形单元的差分方程为 （4.2.55）其中：因子为 （4.2.56）图4.

32、8 矩形差分网格中的井点单元然而，当井孔附近的径向流扩展到相邻节点14时，这些相邻节点上的水头还可以通过式（4.2.54）获取 （4.2.57）式中：为节点i与节点0的距离（图4.8）。把式（4.2.57）带入式（4.2.55）得 （4.2.58）其中： （4.2.59） （4.2.60）其中：系数A、B、E都只取决于网格的形状，可称为几何因子。式（4.2.28）建立了井孔流量、计算井点水头H0和井孔实际水头之间的定量关系，即 （4.2.61）这是综合考虑面源补给、地下水储存量变化和非正方形网格的条件下推导出来的校正公式，因此更具有普遍意义。下面讨论因子的影响。如果井点单元的相邻单元都具有相同

33、的横向尺寸x和纵向尺寸，则有 （4.2.62）因此，式（4.2.61）中关于因子的项被抵消了，校正公式可进一步简化为 （4.2.63）只要在模型网格剖分时稍加注意，这一点很容易实现。实际上，在相邻单元尺寸与井点单元尺寸相差不显著的情况下，几何因子E接近于井点单元的面积，从而使因子对井孔校正结果的影响比较小，基本可以忽略。由于具有这种能够被抵消的性质，即使因子是非均匀分布的，其影响也将在一定程度上削弱。当然，如果因子的分布在井孔附近剧烈变化，假设因子与径距无关还是可能带来一定的误差，有待改进。当井孔附近的网格尺寸完全均匀时，有 （4.2.64）则校正公式可还原为 （4.2.65）其中等效半径与式

34、（4.2.48）相同。这说明Prickett（1967）等提出的校正公式是更普遍的校正公式（4.2.61）的一个特例。有限差分模型中井孔的校正需要区别考虑已知流量井孔和已知水头井孔。对于已知流量井孔，先把井孔流量Qw带入常规差分方程（4.2.55），得到井点单元的计算水头H0，然后调用校正公式（4.2.61）求出实际的井孔水头HW。而对于已知水头的井孔，则先把校正公式（4.2.61）表示的井孔流量Qw带入常规差分方程（4.2.55），其中Hw取已知的井孔水头，利用新的差分方程计算出经典单元的特征水头H0，然后再次调用校正公式（4.2.61）计算出井孔流量Qw。这样就解决了常规有限差分模型的“算

35、不准”问题。关于多边形网格有限差分模型中井孔的校正，将在有限元模型的井孔校正方法中给出。六、差分方程组的求解方法设有限差分模型中又M个节点，对于任一节点p都可以建立形如式（4.2.26）或式（4.2.35）和式（4.2.36）的差分方程，不妨用下述形式的通式描述 （4.2.66）式中：下标p-i（注意中间并非减号，为两者连接号）为与节点p相邻的第i个节点在整体网格中的序号；为节点p和节点p-i之间的水流关联系数；为对应节点p、取决于前一时刻水头分布和本时步之间步长的已知变量。考虑与待求变量无关，则上述方程组为线性方程组，可用矩阵和向量表示为 （4.2.67）式中：为渗透矩阵；为节点水头列向量；

36、为右端项列向量。改方程组的解为 （4.2.68）运用线性方程组的求解方法可得到上述解。对于节点数M较小的问题，线性方程组的求解可采用Gauss消去法，LU分解法等直接法。然而，当节点数很大时，直接法的工作量太大，可考虑采用其他方法，如迭代法、共轭梯度法等。为了说明迭代法的思路，我们把式（4.2.66）中表示时刻的上标k暂时省略，方程组变为 （4.2.69）关于节点p的方程可转化为 （4.2.70）上式提供了一种迭代求解的途径。迭代法的步骤是首先令节点的水头等于前一时刻的水头，即 （4.2.71）式中：上标0为第0次迭代运算。再调用式（4.2.70）进行迭代运算，每次迭代运算的公式为 （4.2.

37、72）式中：上标n为第n次迭代。经过足够多次的迭代运算之后，前后两次迭代得到的水头值几乎相等，判断指标为 （4.2.73）式中：为每次迭代结果与前一次迭代结果的最大差距，如果他小于允许误差，迭代过程即可停止 （4.2.74）其中：是预先设置的收敛标准。为了使迭代收敛更快，可以引入松弛因子，并把迭代公式（4.2.72）改为 （4.2.75）松弛因子的取值范围是12，这种迭代格式称为超松弛法。在采用迭代法求解某些非线性问题时，可能使用才能得到收敛结果，这种迭代格式称为低松弛法。第三节 地下水流有限元模型有限元数值模拟技术在固体力学和流体力学中都有广泛的应用，也经常用于进行地下水流的数值模拟。有限元

38、法的原理与有限差分法有所不同，有限差分法以偏微分方程的差分近似为基础，而有限元法以积分方程的离散近似为基础。在历史上，曾经有两种推导有限元方程的途径：迦辽金（Galerkin）法与里茨（Ritz）法。Galerkin法采用加权余量的积分公式（以平面模型为例）为： （4.3.1）式中：R（H）为近似解H形成的余量函数；w为权函数；为模型的某个子域。Rietze法是利用与地下水偏微分方程等价的函数极限值函数进行求解，即 （4.3.2）式中：是等价范涵，积分内公式F定义为使取极小值的函数H（x，y，t）恰好满足地下水流偏微分方程。于是，地下水方程的近似解通过以下极值条件方程 （4.3.3）得到，其中

39、Hp是模型子域上的节点水头。已有证明，当有限元网格单元采用相同的插值函数及形函数时，上述两种方法实际上是完全等价的，所形成的代数方程组相同。为此，本书中不再对有限元法的积分方程进行推导。一、单元与形函数最常见的平面有限单元网格是三角形网格，而四边形网格也有较多的应用。三位有限元 模型一般采用分层网格建立6节点或8节点等参单元。有限元的积分方程就是在网格单元的基础上近似求解的。图4.9三角形单元及其顶点编号首先讨论平面三角形单元（图4.9）。单元e的三个顶点编号按照逆时针顺序分别为i、j、k，水头在单元内任意坐标点（x，y）的分布特征采用线性插值函数进行近似描述 （4.3.4）式中：1、2、3为

40、几何系数。把节点i、j、k处的坐标和水头代人（4.3.3）得 （4.3.5） （4.3.6） （4.3.7）这是线性方程组，求解得 （4.3.8） （4.3.9） （4.3.10）其中： （4.3.11） （4.3.12） （4.3.13） （4.3.14）式中：为三角形单元面积。将式（4.3.8）至式（4.3.10）代人式（4.3.4），则插值函数可改写为 （4.3.15）其中: （4.3.16）是节点i、j、k在单元e内的形参数。对于固定坐标点（x，y）形函数具有以下的性质 （4.3.17）式中：上标e表示求和只能在同一三角形单元内进行。当所有节点的水头已知时，就可以把形参函数代人式（4.

41、3.15）计算三角形单元内任意坐标位置的水头。图4.10任意四边形单元的坐标变换四边形单元也可以采用类似的方法建立插值函数和形函数，但是在形式上复杂一些。首先，任意的四边形都可以影射为一个正方形，图4.10（a）中（x，y）坐标系下的四边形1-2-3-4，可通过变换，形成4.10（b）中坐标系下的正方形1-2-3-4。这个正方形单元被称为等参单元，其形函数定义为 （4.3.18）式中：为节点i的影射坐标。等参单元内的任意一点，通过上述形参数与实际四边形单元内的对应坐标点建立如下关系，即 （4.3.19）式中：xi、yi为节点i的实际坐标。因此，等参单元内形函数与实际单元函数具有一一对应的关系

42、（4.3.20）而等参单元内的任意一点的水头，也可以通过形函数 变换为坐标下处的水头 （4.3.21）在建立有限元方程时，需要用到一下类型的偏导数 （4.3.22） （4.3.23）， （4.3.24）二、基于三角形单元的地下水流有限方程对于三角形网格中的一个节点p（图4.11），建立有限元方程的积分区域包括了与p点相邻的所有节点p-1、p-2等，即以p点为共同顶点的所有三角形单元所覆盖的区域。积分采用分片发，即在相邻的三角形单元内分别积分，然后求总和，得到对应于p点的有限元方程。图4.11三角形网格中的节点及其水均衡控制区（阴影部分）通过建立有限元网格节点的水均衡控制区，采用水量均衡原理，也

43、可推导出地下水流的有限元方程，与积分法得到的结果相同。这种推导方法更加清楚的反映了地下水流的物理机制。如图4.11所示，节点p的谁均衡控制区是一系列三角形重心与侧边中点的连线所围成的区域，每个三角形单元中，属于节点p谁均衡控制区的面积是其单元面积的1/3。因此，节点p谁均衡控制区的面积为 （4.3.25）式中：e-n为第n个相邻三角形单元的整体编号；Ae-n为其面积。地下水通过控制区周边流入控制区内侧向流量QH，可以用流速在控制区边界线上的线积分得到 （4.3.26）式中：为流速向量；为边界线上微分线段的内法线矢量；和分别为x方向和y方向的流速分量；m为含水层的厚度。设渗透系数的主轴与坐标轴一

44、致，则流速向量可根据达西定律获取， （4.3.27）代人式（4.3.26）有 （4.3.28）确定上述积分的方法是在节点p相邻的三角形单元内分别计算，然后求和，即 （4.3.29）其中，每个三角形单元e-n流量贡献为（4.3.30）式中：S-i为在三角形单元e-i内的积分路径。以图4.11中的单元e-1为例，节点p控制区在该单元内的积分路径为A-B-C，其中A为点p和点p-1连线的中点，B为点p和点p-2连线的中点，C为单元e-1的重心。根据三角形单元内水头的插值函数式（4.3.15），有 （4.3.31）根据式（4.3.16）得 （4.3.32）同理有 （4.3.33）可见水力梯度在三角形单

45、元内是一常量，可以提取到式（4.3.30）的积分号外部，即（4.3.34）这样在线积分中只剩下积分路径在x轴和y轴上的投影，有 （4.3.35）根据点A和点与单元e-1的三个顶点之间的关系，有 （4.3.36） （4.3.37）将式（4.3.36）和式（4.3.37）代人式（4.3.35），得 （4.3.38）将式（4.3.38）代人时（4.3.34），得=（4.3.39）三角形单元e-1中的节点p、节点p-1以及节点p-2分别与图4.9中普遍单元的顶点i、j、k相互对应，有了式（4.3.39）之后，引入如下的系数来反映单元e内节点与节点之间的关系 （4.3.40）式中：Cpl为节点L与节点p

46、的控制区之间得到关系。利用式（4.3.40）可以把式（4.3.39）改写为 （4.3.41）代人式（4.3.39）得到流向节点p控制区的总侧向流量 （4.3.42）式中：p-j和pk分别表示沿着单元e-i从节点p逆时针移动遇到的第一个节点与第二个节点。以图4.11中的节点p为例，展开式（4.3.42），得到（4.3.43）在整个网格中，节点之间在各个单元中的关联系数可以合并为 （4.3.44）式中：Gpl为节点p与节点L之间在整体网络上的关联系数；Epl为共同拥有节点p和L的单元的数目。，说明节点p和L之间没有直接关系。关联系数的重要特征是对称性：Epl=Gpl。运用式（4.3.44），流向节

47、点p控制区的侧向流量可表示为 （4.3.45）式中：N为网格中的总节点数。对于承压含水层，在时间步长#tk内，节点p控制区内水体积的增加量为 （4.3.46）式中：k为时间步长。标准的有限元法是把差值函数引入上述积分，并分别在相邻单元内积分再求和，导出#Vw中包含全部相邻节点的水头变化信息。目前，地下水流有限元模型中，普遍采用储量集中处理法，即用节点p的水头变化表示控制区的平均水头变化，有 （4.3.47）其中：Ap按照式（4.3.25）计算，这种简单的做法反而可以改善有限元的计算结果。根据书均衡原理，控制区水体积的增量与外部流量补给的水量相等，即 （4.3.48）式中：是其他外部不计流量。利

48、用式（4.3.45）和式（4.3.47），可以展开式（4.3.48）为 （4.3.49）这就是承压含水层模型中对应节点p的隐式有限元方程。潜水含水层的饮食有限元方程课写为 （4.3.50）式中：为给水度。关联系数中的导水系数Tx、Ty与当前时刻的带球水头有关，因此属于非线性方程。三、三角形四边形混合单元地下水流有限元方程对于三角形单元，已经通过式（4.3.40）得到了单元内e内关联系数Gpl的表达式。在有限元模型中，对于任意单元体，节点之间关联系数的积分形式为 （4.3.51）这个积分表达式对三角形单元同样是有效的。对于四边形单元，需要将上述积分转换到射影的等单元上进行，利用以下性质 （4.3

49、.52）式中：为坐标变换的Jacob矩阵，而为其行列式，有 = （4.3.53）并且 （4.3.54） = （4.3.55） 而x、y对、的偏导数具有式（4.3.22）的形式。利用式（4.3.52），积分式（4.3.51）可转变为（4.3.56）注意积分符号内的偏导数和也是影射坐标（，）的函数。 上述积分可采用高斯积分公式进行计算（4.3.57）其中：（i，j），i=1，2，3，4为四个高斯积分点，并有 （4.3.58） 在确定了单元内的节点关联系数之后，就可以利用（4.3.44）组装整体关联系数。而出水量部分仍然采用储量集中处理法，节点p控制区在单元e内所占的面积为 （4.3.59） 而节点p控制区的总面积是（4.3.60）其中：为与p相邻的单元数。 最后可以得到承压水层三角形四边形混合单元模型的有限元方程为（4.3.61）四、地下水流有限元模型的孔井校正与有限差分模型一样，地下水流有限元模型在遇到井孔时也需要进行校正。井孔的校正思路与有限差分模型中所采取的思路是相同的即引入井孔附近的径向流解析解在研究地下水流有限差分模型的