第5章流网

校正：解千里，骆慧斌

统稿：吴小龙

5.1流网图构建

我们已在第2章看到，地下水流系统可以用三维等势面集和一个对应的正交流线集来表示。如果通过三维系统能够挑选出一个有代表性的二维剖面，那么等势线和流线的集合就形成了“流网”。构建流网图是地下水流分析中最有力的工具之一。

在2.11节和图2.25中，我们已看到，流网可被当作一个有界二维稳定流的解。这个解，需要知道地下水流经的区域，这个区域边界上的边界条件，和这个区域内水力传导率的空间分布。在附录III中提供了对应的数学解析解法。在这一节中，我们将学习如何使用图像方法构建流网，而不使用复杂的数学方法。

均质、各向同性系统

让我们首先考虑一个均质、各向同性、并充分饱和的水流区域。对于这样一个区域中的稳定流来说，可能存在三种型式的边界：(1)不透水边界，(2)定水头边界，和(3)地下水位边界。我们首先考虑一个不透水边界附近的水流[图5.1(a)]。由于这里没有水流能够穿过边界，所以靠近边界的流线必然是平行于边界的。同时，等势线必然是以直角与边界相交。引用达西定律，并使穿过边界的单位面积流量等于零。这样，我们就得到了边界条件的数学解释。对于在xz平面上平行于轴的边界来说：

$\frac{\partial h}{\partial x} = 0$ 或 $\frac{\partial h}{\partial z} = 0$ (5.1)

**图 5.1** 在（a）不渗透边界，（b）定水头边界，（c）地下水位边界周围的地下水流。

实际上，一个流网中的任何流线都构成一个假想的不透水边界，即没有水流能够“穿过”流线。在构建流网图时，往往需要缩小水流区域的尺寸，方法是只考虑对称部分的一侧或另一侧。如果明确对称线也是一条流线，那么对称线可以被视为边界，其边界条件就如方程(5.1)所示。

水力水头为常数的边界就是一条等势线[图5.1(b)]。流线必须以直角与边界相遇，并且与其相交的等势线必须与边界平行。其数学表达是

h = c (5.2)

在地下水位上，其压力水头ψ等于零，同时边界条件根据水头关系式h = ψ + z变为

h = z (5.3)

正如图5.1(c)所示，对于一个地下水补给情况来说，地下水面既不是一条流线，也不是一条等势线。它只是一条已知h但高度变化的线。

如果我们知道一个均质、各向同性的水流区域内地质材料的水力传导率K值，就可能用流网计算出通过这个系统的流量。图5.2是曾在图2.25(a)中提过的一个简单情况的完整流网。两条相邻流线之间的面积称为“溪管”或“流管”。如果各流线是等距的，则通过每个流管的流量是相同的。考虑图5.2中通过区域的水流。如果距离AB和BC分别是ds和dm，AD和BC之间水力水头下降值是dh，则穿过此区每单位深度(垂直于纸面方向)横剖面的流量是

$dQ = K\frac{dh}{ds}dm$ (5.4)

在稳定流的条件下，通过流管内单位深度的任何剖面(指的是AD, EH, 或FG)而流过的流量必然也都是dQ。换句话说，通过一个流管任何部分的流量，可以根据流管的一个部分计算得出。

如果我们任意确定一个方形流网，使其ds = dm，则方程(5.4)变为

$dQ = K dh$ (5.5)

对于一个具有m个流管的水流系统来说，其总的流量是

$Q = mK dh$ (5.6)

如果穿过水流区的总水头下降值是H，同时在流网中有n个水头分段(H = n dh)则

$Q = \frac{mKH}{n}$ (5.7)

对于图5.2来说，m = 3, n = 6, H = 60米，同时从方程(5.7)可知，Q = 30K。对于K = 10-4米/秒来说，Q = 3 × 10^-3 m³立方米/秒(垂直于流网方向上的剖面)。

使用方程(5.7)时必须谨慎。它只能用于具有一个补给边界和一个排泄边界的简单水流系统。至于更复杂的系统，最好是先计算一个流管的dQ，再以流管数目乘之来获得Q值。

图5.3所示的流网包含通过坝基地基的渗流，地基之下是一个不透水边界。我们可以根据这张流网表明三个论点。

即使是最简单流网中的“方格”，实际上也都是曲线方格；在这些方格的中心，数值是相同的。换句话说，它们实质上是一个圆圈，这个圆圈对于其全部四个边界线来说都是正切的关系。
并不需要流网的全部边界都是有限边界；像图5.3中水平无限伸展的地层那样，当水流区域在一个或更多方向上无限伸展，也是可以处理的。
—个流网，可以在边缘上用一个“局部的”流管来处理。

对于图5.3中的流网来说， $m = 3\frac{1}{2}$ 。如果H = 100米，K = 10^-4米/秒，那么，由于n = 6,我们得Q = 5.8×10^-3立方米/秒(垂直于流网的剖面)。

在均质、各向同性的介质中，水力水头的分布仅仅决定于边界条件的形态。流网的定性与介质的水力传导率无关。仅仅当进行定量的流量计算时，水力传导率才有用。也应该注意到流网是无量纲的。图5.2和图5.3中的流网，不管你把它看作是几米大小的方格，还是上千米大小的方格，它们都是同样有效的。

勾划流网是一种艺术性的工作，人们经常是在用试错法多次反复的基础上进行的。一些水文学家已极为熟练，能够很快使流网达到可采纳的程度。对于其他人来说，它还是个不断经受挫折的起点。对于一个均质、各向同性介质中的流网来说，曲线构图的规则是十分简单的。我们可把它总结为如下几点：(1)在整个系统内，流线和等势线必须以直角相交；(2)等势线必须以直角与不透水边界相连；(3)等势线必须平行于不变水头边界；(4)如果流网所划出的方格是在野外某个局部作出的，那么，除去在边缘处可能有部分流管特殊外，这种方格必须对整个野外测区具有代表性。

非均质系统和正切定律

当地下水流线穿过具有不同水力传导率的两个层组之间的地质边界时，它会像光线通过一种介质进入另一介质时的表现一样发生折射。但是，与斯涅尔定律(是一个正弦定律) 有区别的是，地下水流线的折射遵守正切定律。

考虑图5.4中所示的流束，水流从水力传导率为K₁的介质向水力传导率为K₂的介质流动。其中流线具有与书页垂直的单位深度，其角度和距离已在图上指出。对于稳定流来说，入流量Q₁，必然与出流量Q₂相等。或者，根据达西定律

$K_1a\frac{dh_1}{dl_1} = K_2c\frac{dh_2}{dl_2}$ (5.8)

其中dh₁是流过距离dh₁时的水头下降值，dh₂是流过距离dh₂时的水头下降值。在dh₁和dh₂中都以相同的两个等势线为界，因此dh₁=dh₂；同时，根据几何关系，a = b cos θ₁且c = b cos θ₂。注意到b/dl₁ = 1/sinθ₁且 b/dl₂ = 1/sinθ₂方程(5.8)变为

$K_1\frac{\cos \theta_1}{\sin \theta_1} = K_2\frac{\cos \theta_2}{\sin \theta_2}$ (5.9)

或

$\frac{K_1}{K_2} = \frac{\tan \theta_1}{\tan \theta_2}$ (5.10)

方程(5.10)构成了非均质介质中地质边界地下水流线的正切折射定律。已知K₁, K₂和θ₁，人们就可以根据方裎 (5.10) 求θ₂。图5.5表示具有K₁/K₂ = 10的两种地层间的流线折射。根据它们本身的性质，流线将倾向于以高透水性层组为水流通路，而试图以最短的路径通过低透水性层组。因此，在含水层-滞水层 (在透水性上相差2个数量级或更多) 系统内，含水层中的流线将变为几乎是水平的，而在滞水层内，流线则几乎是垂向的。当人们考虑到列于表2.2中的水力传导率的广大范围时，就会很清楚地知道，相差2个数量级或更多是十分常见的情况。

**图 5.5** 流线在分层系统中的折射现象（根据Hubbert, 1940）。

如果人们试图在图5.5的曲线上描绘出等势线以完成水流系统，将会很快明白，在所有层组中构成方格是不可能的。在非均质系统中，一个层组中的方格在其它层组中变成了长方形。

我们可以把非均质、各向同性系统中曲线流网构图的规则总结如下：(1)流线和等势线必须在整个系统中以直角相交；(2)等势线必须以直角与不透水边界相连；(3)等势线必须与不变水头边界平行；(4)在地质边界上必须满足正切定律；(5)如果流网所画的方格是在一个层组中的一部分中画出的，那么这种方格必须是在整个层组中都存在。同时，在所有水力传导率相同的层组中也都存在。在水力传导率不同的层组中所画出的将是长方形。

最后两条规则使在复杂的非均质系统中极难画出精确的定量流网。但是，定性的流网，只保持其中的正交性而不尝试构建方格，对于了解地下水流系统仍是一个巨大的帮助。图5.6是原已在图5.3中介绍过的在一个坝下渗流问题中的定性示意流网。现在增加了一层基岩。

各向异性系统和变换剖面

在均质但却各向异性的介质中，流网中流线与等势线不呈直角而显得更加复杂化。Maasland(1957)，Bear和Dagan(1965)，Liakopoulos(1965b)对这一现象中的理论原理进行了研究。同时Bear(1972)提出了丰富的理论观点。在这一节中，我们将主要检验之前较少被考虑的非正交性条件下的实际应用。这将会涉及使用“变换剖面”对流网进行构建。

例如，我们考虑一个均质、各向异性二维区域中的水流。若主轴方向上的水力传导率分别为K_x和K_z，则水力传导率椭圆(图5.7)将有半轴 $\sqrt{K_x}$ 和 $\sqrt{K_z}$ 。让我们变换一下水流区域的比例，使新的坐标系统X和Z以下式与xz系统相关联

$X = x$
(5.11)

$Z = \frac{z\sqrt{K_x}}{\sqrt{K_z}}$

若K_x>K_z，这个变换将使水流区域在垂向方向上规模扩大，也使水力传导率椭圆扩大为一个具有半径为 $\sqrt{K_x}$ 的圆(即图5.7中的外圆)。同时，这个虚构的扩大了的水流区域将起到和水力传导率为K_x的均质系统中一样的作用。

**图 5.7** 各向异性介质（*K_x*/*K_z*5）的水力传导椭圆。两个圆表示两种可能的各向同性变换。

这个变换也可在稳定流方程的基础上保持有效性。例如在原来xz座标系统中，根据方程(2.69)，对于一个各向异性介质来说，我们有

$\frac{\partial}{\partial x}\left(K_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial z}\left(K_z\frac{\partial h}{\partial x}\right) = 0$
(5.12)

两边都以除之K_x得

$\frac{\partial^2 h}{\partial x^2} + \frac{\partial}{\partial z} \left(\frac{K_z}{K_x}\frac{\partial h}{\partial x}\right) = 0$

根据方程(5.11)中的第二个表达式，对于变换了的剖面来说，我们得

$\frac{\partial}{\partial z} + \frac{\sqrt{K_x}}{\sqrt{K_z}}\frac{\partial}{\partial Z}$
(5.14)

注意到方程(5.11)中的第一个表达式，将方程(5.14)的运算应用到方程(5.13)的两项微分式中，我们得

$\frac{\partial^2 h}{\partial X^2} + \frac{\partial^2 h}{\partial Z^2} = 0$ (5.15)

这是在变换剖面中，对于一个均质、各向同性介质的水流方程。

同样的，我们也可以按照如下关系进行变换

$X = \frac{x\sqrt{K_z}}{\sqrt{K_x}}$
(5.16)

$Z = z$

水流区域则会在x方向上缩小。在这种情况下，水力传导率椭圆将被变换为图5.7中的小圆。同时，这个虚构的水流区域将起到和水力传导率为K_z的均质系统中一样的作用。

在变换剖面概念的应用中，均质、各向异性流网的构图步骤是不证自明的：(1)用方程(5.11)或(5.16)进行座标变换；(2)在假想的、已变换的剖面中按照均质、各向同性介质的规则构出一个流网；(3)逆变换坐标比率。

**图 5.8** （a）均质、各向异性 $\sqrt{K_x}/\sqrt{K_z} = 4$ （b）在变换的各向同性剖面中的流网，（c）在原始各向异性剖面中的实际流网，T——变换；I——逆变换。

**图 5.8** （a）均质、各向异性 $\sqrt{K_x}/\sqrt{K_z} = 4$ （b）在变换的各向同性剖面中的流网，（c）在原始各向异性剖面中的实际流网，T——变换；I——逆变换。

图5.8是这项技术的一个实例。图5.8(a)中所示的边界值问题代表从h = 100 米处的地面水池向h = 0米处排水水流的一个垂向剖面。这种排水水流代表了在同一深度垂直纸面的很多平行的排水流网之一。因为整个水流系统存在对称性，垂向的边界可以被视为“假象的”不透水边界。下部边界则是一个实际边界，它代表土壤的基底。其下面的土质或岩石层组的水力传导率比其上面的低若干个数量级。如果将纵轴z = 0设在排水处，地面水池设为 z = 100，那么，根据h = ψ + z和已给的h值，我们在两个边界处都得到ψ = 0。在地面水池处，这个条件意味着土壤是饱和的，水池是起点，它的深度是0。在排水点，ψ = 0意味着自由流动的条件。水流场中的土壤有一个K_x/K_z的各向异性水力传导率。所以图5.8(b)的变换剖面有一个 $\sqrt{K_x}/\sqrt{K_z} = 4$ 的垂向扩大。图5.8(c)表示逆变换的结果，其中由变换剖面作出的均质、各向同性流网又重新变回为实际尺寸的水流区域。在逆变换时，图5.8(b)中任何一点(X, Z)的水力水头都变为了图5.8(c)中点(x, z)处的水力水头。

很明显，变换剖面的尺寸取决变换时所用的方程(5.11)或(5.16)。但是在任何一种情况中，流动区域的形状，和构成的流网，都是一样的。

如果所求的是排泄量或水流速度，那么通常最容易的方法是在变换剖面中进行这些计算。因为这样，问题就变为了单纯使用水力传导率计算的问题了。很清楚，对于一个垂向扩大了的剖面使用K_x值，或是对于一个水平缩小的剖面使用K_z值都是不正确的。根据从图5.7所作出的推断，对于这种同一问题的两个等同的表达式来说，将产生两组不同的定量计算成果。事实上，所用的正确值应是

$K' = \sqrt{K_x \cdot K_z}$ (5.17)

方程(5.17)的有效性在于：可以使得两个等同的变换所产生的水流区域都是相等的。要证明它，需要把达西定律应用在每个变换剖面中的单一流管中。

各向异性在地下水流网性质上的影响，可以于图5.9说明。它和之前提到的图5.8是相同的边界问题。各向异性流网[图5.9(a)和5.9(c)]的一个最重要特征是它们缺乏正交性。本节中所引出的变换技术对这一现象提供了非直接但却还能满意的解释。

**图 5.9** 中水流问题的流网，其中 $\sqrt{K_x}/\sqrt{K_z} =$ 分别等于（a）1/4，（b）1，（c）4（根据Masland, 1957)。

**图 5.9** 中水流问题的流网，其中 $\sqrt{K_x}/\sqrt{K_z} =$ 分别等于（a）1/4，（b）1，（c）4（根据Masland, 1957)。

人们可能在很多情况之下都会需要作流网图(以野外地下水位数据为基础)。如果地质层组是已知的各向异性，那么在依据等势数据去推断水流方向时，必须十分谨慎。如果需要一个完善的流网，自然就需要用变换剖面去完成。但是如果仅仅需要各个特定点处的水流方向，那么就还有另外一种构图方法可用。如图5.10虚线代表xz野外场中所感兴趣的点的等势线的方向趋势。围绕这一点构出一个“反转的”水力传导率椭圆。这个椭圆的半轴为 $1/\sqrt{K_x}$ 和 $1/\sqrt{K_z}$ (而不是图5.7中的 $\sqrt{K_x}$ 和 $\sqrt{K_z}$ )。若水力梯度的方向线和椭圆相交于点A，在A点画出一条相对于椭圆的正切线，那么水流方向就是垂直于这条正切线的一条线。作为这种构图法的一个实例，人们可以把图5.10中的成果与图5.9(c)比较。

**图 5.10** 在各向异性 $\sqrt{K_x}/\sqrt{K_z} = 5$ 围4个值的平均数区域中水流方向的确定。

**图 5.10** 在各向异性 $\sqrt{K_x}/\sqrt{K_z} = 5$ 围4个值的平均数区域中水流方向的确定。

5.2类比模拟流网

对于一个xz座标系统中均质、各向同性介质中的水流来说，一个流网中的等势线代表了的该区域稳定流的边界值问题的解h(x, z)。流网的便是拉普拉斯方程的一个间接解

$\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial z^2} = 0$ (5.18)

这个方程是数学物理中最常见的偏微分方程之一。可以用于描述穿过固体的热流和穿过导体介质的电流等其它物理现象。对于电流的情况来说，拉普拉斯方程采用如下形式

$\frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial z^2} = 0$ (5.19)

其中V是电势或电压。

方程(5.18)和(5.19)的相似性，显示了电流与地下水流之间在数学与物理模拟上的相似性。这两个方程都是根据线性流定律在不同的情况下推导出的。它们都基于一种连续性关系，即对于达西定律而言流体质量守恒，对于欧姆定律而言电荷守恒，试比较一下欧姆定律

$l_x = -\sigma\frac{\partial V}{\partial x}$ (5.20)

和达西定律,

$v_x = -K\frac{\partial h}{\partial x}$ (5.21)

我们就发现它们的相似性。单位面积流量v_x(即每单位面积的排泄量)与电流密度I_x(即每单位面积的电流)相似；水力传导率K与单位电导性σ相似；水力水头h与电势V相似。

电流与地下水流之间的相似是两种类型模拟模式的基础(这些模拟模式已被证实对于绘制定量流网是有用的)。第一种类型利用电导纸，第二种类型使用电阻网。

电导纸模拟

让我们再次考虑初次表示在图5.8，现在又重新出现在图5.11(a)中的水力问题。电模拟拟[图5.11(b)]由一张剪成地下水流场形状的电导纸构成。电力的供应是通过在边界造成一个电位差，同时使用伏特计的探针来测量整个电导纸上的电势分布的。图5.11(b)上所示的V = 100的定水头边界，是用高电导性的银涂料作成的。不透水边界是使用边缘不相互接触的纸模型来模拟的。这种方法通常能有效的找出等势线的确切位置。因此能够很快地完成一个完善的等势网。

这种方法仅适用于二维均质、各向同性系统中，但它却能解决复杂的区域形状和边界条件。商店中现成电导纸在电导性上存在差别，会导致随机误差，限制了这个方法在定量上的精确性。这个方法的两次最细致的应用是Childs(1943)在近地表排水系统的理论分析和Toths(1968)在阿尔伯塔对一个野外地区的区域性地下水流网的研究。

**图 5.11** 电类比模拟流网：（a)稳态水文地质边界问题；（b)电导纸类比；（c)电阻网类比。

电阻网模拟

在电模拟中利用电阻网，是以电导纸模拟的同样原理为基础的。在这个方法中[图5.11(c)]，水流场被一个在网格节点上相互联结的电阻器网所代替。电通过每个电阻器的流动与地下水通过一个流管的水流相似(这个流管平行于电阻器，并有一个由电阻器间距乘以单位深度反映出来的横剖面面积)。对于通过一个单独电阻器的电流来说，方程(5.20)中的I必须被看作是电流，并且σ是等于 1/R的，其中R是电阻器的电阻。像在电导纸模拟中一样，穿过模型的定水头边界设立一个电位差。用一个探针去确定网中每个节点的电压。当把这些值记录下来并且连成等值线时，等势网就画出来了。

改变网中的电阻，就能够用电阻网模拟方法去分析非均质、各向异性的系统了。它的精确性和多用性要优于电导纸模拟。但是它没有数值方法(将在下节介绍)那样的灵活性。

Karplus(1958)提供有类比模拟的详细手册，电阻网模拟曾被Luthin(1953)在排水工程中用来描画地下水流网，被Bouwer和Little(1959)应用于饱和-非饱和系统中。Bouwer(1962)曾用这个方法去分析在补给水池下面形成的地下水圆丘的形态。

在地下水文学中电模拟方法应用最广泛的地方，是在含水层中非稳定水流分析中所使用的电阻-电容网方法。这种应用将在8.9节中讨论。

5.3数值模拟流网

构成一个流网的水力水头场h(x, z)能够通过两种途径从有关稳定流边界值问题中用数学推导出来。第一种途径是根据2.11节和附录III中所叙述的解析方法，第二种途径是数值方法。解析方法适用于在水流区域、边界条件和地质形态比较简单与规则的水流问题中，而我们将在本节中所看到的数值方法则是更全能的，但是它的应用毫无疑问地是要与数字计算机的应用结合在一起。

数值方法是近似方法。它是以构成水流区域的连续体被离散化为基础的。在离散化中，区域被分割为有限数量的块体，每个块体都具有它自己的水文地质特性。每个块体中心有一个节点，该节点处的水力水头是由整个块体来确定的。图5.12(a)表示一个长方形水流区域的7 × 5个块体的中心节点网格(在x方向上从i = 1到i = 7,在z方向上从j = 1到j = 5)。

现在让我们来考察内部节点之一附近的水流状态——比如节点块体i = 4, j = 3处及其周围相邻块体的水流情况。为了简化我们以图5.12(b)所示重新对节点标注。如果水流是从节点1到节点5发生的，我们就能够根据达西定律去计算垂直于纸面单位深度横剖面的流量Q₁₅：

$Q_{15} = K_{15}\frac{h_1 - h_5}{\Delta z}\Delta x$ (5.22)

当假设在每一种情况下水流都是直接流向中心节点时，我们就能够写出对于<Q₂₅, Q₃₅和Q₄₅的类似表达式。对于稳定流来说，按照流体质量守恒定律，这四股水流之和必须为零。如果介质是均质和各向同性的，则K₁₅ = K₂₅ = K₃₅ = K₄₅，如果我们任意选择一个正方形的节点网格，以使Δx = Δz，这四项之和将是

$h_5 = \frac{1}{4}(h_1 + h_2 + h_3 + h_4)$ (5.23)

这个方程就称为“有限差分方程”，如果我们再转化为图5.12(a)的ij标识，方程5.23就变为

$h_{i,j} = \frac{1}{4}(h_{i,j-1} + h_{i+1,j} + h_{i,j+1} + h_{i-1,j})$ (5.24)

**图 5.12** (a)数值流网模拟中的块体中心节点网络 (b)内部节点及其相邻的节点块体 (c)对内部节点、基底不渗透边界节点、以及不渗透角端节点的有限差分星式示意图。

这个形式中，方程(5.24)节点网格中所有内部节点都成立。它优雅地说明了一条朴素的真理：在稳定流中，在均质、各向同性介质中，其任一节点处的水力水头是其周围4个值的平均数。

对于沿基底边界(并假设是不透水边界)的上一个节点的有限差分方程，我们可同样写出如下形式：

$h_{i,j} = \frac{1}{4}(h_{i-1,j} + h_{i+1,j} + 2h_{i,j+1})$ (5.25)

以及对于不渗透角端节点的

$h_{i,j} = \frac{1}{4}(2h_{i-1,j} + 2h_{i,j+1})$ (5.26)

方程(5.24)到方程(5.26)的图解示意在图5.12(c)的三个有限差分星式示意图中。

简而言之，在节点网格中，对每一个节点建立一个有限差分方程都是可能的。如果有N个节点，就有N个有限差分方程。也有N个未知数—在N个节点上的N个h值。因此我们在N个未知数中产生了N个线性代数方程。如果N值不大，我们就能够利用像Cramer法则那样的技术直接解出方程。但是如果N很大(一般在数值流网模拟中N都很大)，我们就必须引入一种称为“逐次渐近法”(或“松弛法”)的方法才更有效率。

让我们继续思考图5.11(a)所示的水流问题，并假设我们要用数值方法建立它的流网。在图5.12(a)的节点网格中，水头已知的节点都划了圆圈：在i = 1, j = 3处h = 0，在j = 5那一排上所有节点的h = 100。松弛法通过对整个节点网多次重复地扫描：从顶到底，从左到右(或者采用任何适当的方式)，在每个水力水头未知的节点上应用相应的有限差分方程。该方法必须在每个节点处假设一些起始的h值。对于当前问题来说，我们可以假设h = 50作为所有未知节点上的起始值，每次在有限差分方程的应用后把最新算得的h值应用在节点上。每将整个系统的节点计算一次，就称为一次“迭代”。在每次迭代之后，所算得的h值就更接近其答案一步。在任一节点，两次连续迭代之间h值的差值称为“残差”。系统中的最大残差将随着迭代的进行而逐渐减少。当最大残差减小到低于预先确定的“允许公差”以后，算法就获得了问题的一个(近似)解。

读者如果需要测验对松弛的理解程度，可在流网的左上部进行一些迭代。如果在节点i = 2, j = 4处把初始值确定为50，那么在第一次迭代之后的h值将是62.5，第二次迭代之后为64.0，在若干次迭代之后所达到的最终值在65和66之间。

数值模拟是能够处理任何形状的水流区域和任何分布状况的边界条件的，它能轻易地对长方形(Δx≠Δz)和非均质、各向异性的水力传导率分布(K_x和K_z和值在每个节点都不同)的节点网格建立有限差分方程。在方程(5.22)中，常用的平均方法将使K₁₅ = (K₁ + K₅)/2(其中K₁和K₅在这种情况中被当作是节点1和5处的垂向水力传导率，它们可能彼此不同，也可能与这些节点上的水平水力传导率不同)。数值模拟也能在图解构造或解析解十分复杂的情况下构建出流网。数值模拟几乎永远是编程为数字计算机程序的，并且计算机程序经常是写成普适的形式，以处理广泛的、多种多样的水流问题，而且仅仅需要新的数据。这与电阻网模拟相比显然是一个明显的进步。因为电阻网模拟在进行一项新的模拟时，需要把原来已拼成的硬件全部拆掉。

本节中所提出的有限差分方程的推导是颇不正规的。它完全可以从拉普拉斯方程开始，并且进一步循着数学推导达到同样的结果。在附录VI中，我们沿着这条思路提出了主要的推导步骤。值得注意的是前面用达西定律和连续性方程到有限差分表达式的非正规推导与2.11节中拉普拉斯方程的推导有同样的两个步骤。

根据Shaw和Southwell, 1941命名的松弛法(“逐次渐近法”)有几个广为人知的别名：“高斯-塞德尔方法”，“利布曼方法”，和连续替代法(method of successine displacements)。它在很多解算有限差分方程组的有效方法中，是一种最简单、但远非最有效的方法。例如，如果在逐次渐近期间，对所计算出的h值按照

$h_{corr}^k = \omega h^k + (1 - \omega)h_{corr}^{k-1}$ (5.27)

进行修正(其中h^k是在第k次迭代时所得的计算值， $h_{corr}^{k-1}$ 是前一次迭代中修正后的值)，那么，这就是另一种新的方法，称为“连续超松弛”法。该方法能显著减少达到收敛的近似解的迭代次数。参数ω称为“超松弛参数”，同时它必然在1≤ω≤2的范围内。

有很多参考书可以很好地为数值模拟专家服务。McCracken和Dom(1964)在他们的Fortran手册中，对计算机模拟提供有初步介绍。Forsythe和Wasow(1960)在一个更高的数学水平上提出过他们的观点。Remson(1971)等人特别对有关地下水流的数值技术作了广泛地讨论。Pinder和Gray(1977)对这类问题作了更髙的水平的处理。

数值方法是Stallman(1956)在一次区域性地下水位的分析中引入到地下水文学文献中来的。Payers和Sheldon(1962)是第一批倡导在区域性水文学的研究中引用稳定流数值模拟的学者。Remson(1965)等人把这个方法用于拟建水库对砂层含水层中地下水位影响的预测中。Freeze和Witherspoon(1966)在他们区域性地下水流的理论研究中创制出了很多数值流网。这个方法已经很早就在农业排水区(Luthin和Gaskell, 1950)和土坝渗漏图型的绘制中得到了广泛地应用(Shaw和Southwell, 1941)。

近几年来，有限差分法已和另一个称为“有限单元法”的数值解算方法同样卓有名望。有限单元法也建立N个可由逐次渐近法求解的方裎(有N个未知数)。但是有限单元法中的节点是一个不规则三角形或四边形的顶点，这些不规则的形状是为每种特殊用途的模拟而设计的。它与有限差分法的正规方形不同，在很多情况下需要更少的节点网格，节省更多的计算工作。有限单元法还能处理有限差分法不能处理的情况。有限差分法要求各向异性层组中的各向异性的基本方向平行于座标方向，如果在水流场中有两个各向异性层组，各有不同的基本方向，有限差分法就陷入了困境，而有限单元法则能够求解。有限单元方程的推导需要高深的数学知识，已超过了这本介绍性的书籍范围，有兴趣的读者可去参考Pinder和Gray(1977)的著作。

有限差分法和有限单元法都已广泛用作为地下水含水层中“非稳定流”数字计算机模拟的基础，这些用法在8.8节中讨论。

5.4饱和-非饱和流网

这里有另外一种非常难于用图解的方法来构成的流网型式。对于同时包括饱和与非饱和两种水流的问题来说，一般都是用数值模拟来导出其静态流网。比如图5.13中所示的流网：它与我们在前节中已多次反复分析过的问题是一样的，即在三个边都是不透水边界的系统中有流向排泄点的水流。但不同的是，它在上部边界表示有小于大气压力的压力水头的垂向刻度，这意味着在地表面土壤是非饱和的。然而如果要使排泄点有出流发生，它在深处就必须是饱和的。图5.13中的定性流网是为一种土壤而建立的，它的非饱和特性曲线已表示在插图中。这些作为ψ的函数的水力传导率K和含水量θ曲线，是从图2.13引来的湿化曲线。

**图 5.13** 均质、各向同性土壤中饱和-非饱和流网。插图所示为土壤中的非饱和特性曲线。

正如在图2.12中所示意说明的那样，在一维饱和-非饱和情况下，从数值模拟的静态饱和-非饱和二维流网中有三种类型的输出：第一是水力水头 h(x, z)图形，能够构成等势网(图5.13中的虚线)；第二是压力水头ψ(x, z)图形(图5.13中的点线)，它是确定水位位置(ψ = 0的等压线)的特定值；第三是含水量θ(x, z)图形，它可在土壤 θ(ψ)曲线的帮助下从ψ(x, z)图形确定出来。例如，沿着图5.13中的ψ = -50 厘米的点线，含水量θ是27%。

流线和等势线在整个饱和-非饱和区形成一个连续的网。它们在整个系统中呈直角交叉。在均质、各向同性、饱和部分，用曲线方格可以划出定量流网，但是它们的流管将随着它们在非饱和带中的穿流而不能呈现为方格形，甚至在均质、各向同性的土壤中也是这样。随着压力水头(和含水量)的降低，水力传导率也降低。同时，为了通过给定流管给出同样的流量，就要有一个增大了的水力梯度。这种现象可以在图5.13流网的上部左角的流管中观察到。该处的水力梯度是向地表而增加的。

Luthin和Day(1955)在水文文献中引入了一个完整的饱和-非饱和概念。他们用数值模拟和一个试验砂箱来导出它们的h(x, z)图形。Bawer和Little(1959)用电阻网去分析与图5.13中所示的概念一样的冰碛土排水和地下灌溉问题。由于饱和-非饱和流网适用于溪流发生的问题(第6.5节)，所以在解释悬滞地下水位(图2.15和图6.11)和了解山坡上的水文地质条件时是需要的。Reisenauer(1963)、Jeppson和Nelson(1970)用数值模拟去考察水池和灌渠下面的饱和-非饱和势态。他们的解算法已应用于地下水人工补给的分析中(第8.11节)，Freeze(1971b)研究了饱和带通过土坝渗漏方面的影响问题(第10.2节).

5.5渗透面和裘布衣水流

渗透面、渗出点、自由表面

如果一个饱和-非饱和水流系统存在于一个自由出流边界附近，诸如一条河岸或者一座土坝的下游面，就会在出流边界上出现一个“渗透面”，图5.14(a)中的BC是一个不变水头边界，DC是不透水边界。如果在地表没有水源，AB也将起不透水边界的作用。地下水位EF与出流边界在“渗出点”E相交。所有的水流都必须通过渗出点E下面的渗透面ED而流出该系统。在E点上面，沿AE线，非饱和压力水头ψ是低于大气压的，所以出流到大气圈中来是不可能的，在实际上AE起到了不透水边界作用。ED上的边界条件是 h = z，与地下水位上的条件是一样的。对于这种情况，构建流网的问题在于：在出流边界上分开两种边界条件的渗出点的位置是未知的、预设的。在数值模拟中，对于渗出点的位置需要提供给一个起始的预测值。然后通过一系列反复的稳定流解算的试错演算，才能把渗出点确定下来。

**图 5.14** 自由出流边界渗透面构建。(a) 中饱和-非饱和流网；(b) 自由表面流网；(c) 裘布衣-福尔海默流网。

在饱和-非饱和状态中进行定量流网的构图，需要关于土壤饱和水力传导率K和非饱和特性曲线 K(ψ)方面的知识。在很多工程实践中(包括穿过土坝渗漏的分析)，后一种数据很少有现成的。在这种情况下，经常就去假设通过系统内非饱和带部分的水流是可忽略不计的，或者从另一方面看，当含水量小于饱和度时，土壤的水力传导率与饱和水力传导率对比起来是可忽略不计的。在这种情况下，流网的上部边界就变成了地下水位，并且这个地下水位本身也变成了一条流线。在这种特殊环境下，这个上部边界就称为“自由表面”。以自由表面为边界的饱和系统中的流网可以用通常的方法来构图。但是这里有一个复杂的问题，整个自由表面的位置是未知的推测点，它并不一定恰好在渗出点。在自由表面上的边界条件必须同时满足地下水位(h = z)和流线(等势线必须与它直角相交)两者的要求。其位置经常是用图解试绘与误差的方法来确定。关于工程渗漏问题方面的参考书，诸如Harr(1962)和Cedergren(1967)，提供有图解构图的启示并包括很多稳定、自由表面流网的实例。

图5.14(b)是与图5.14(a)中所示的饱和-非饱和流网等同的自由表面流网，看了这两个图后，我们更确信把地下水位作为流线，对这种特殊的水流系统来说，是一种十分好的近似。出流边界-ED仍然称为渗透面。当我们研究山坡水文学(第6.5节)和考虑土坝下的渗漏时(第10.2节)，我们将在实际中遇到渗透面的问题。

自由表面水流的裘布衣-福尔海默原理

对于非封闭系统中以自由表面为边界的水流问题，经常引入由裘布衣|(1863)所首创并由Forchheimer(1930)所改进的一种方法。这种方法是以下面两点假设为基础的：(1) 假设流线为水平的而且等势线是垂向的；(2) 假设水力梯度是自由表面的坡度，并且不随深度而改变。图5.14(c)中展示了与图5.14(a)同一问题的实质上带有裘布衣假设的等势网，据此我们可以看出，严格的流线构图不再可能了。这种自相矛盾的情况使裘布衣-福尔海默理论被认定为对于实际水流场的一种经验性的近似。实际上，这个理论忽视了垂向水流成分，在实践中，它的价值在于为了分析的目的而把二维系统降低为一维系统。当自由表面的坡度很小而非封闭水流场的深度较浅时，以裘布衣假设为基础的计算比起其他以更严格方法为基础的计算来还是比较有利的。在图5.14(c)中，通过垂直于书面每单位宽度的流量Q由下式给出

$Q = Kh(x)\frac{dh}{dx}$ (5.28)

其中h(x)是在x处水流系统基面以上的自由表面的高度。梯度dh/dx是由x处自由表面的坡度Δh/Δx来给出的。对于稳定流来说，Q在整个系统中都必须是常数，并且这只能在自由表面是一条抛物线时才是真实的。

在均质、各向同性介质中，关于裘布衣-福尔海默理论水流的方程，可以从连续性关系式dQ/dx = 0来进行推导。根据方程(5.28)得

$\frac{d^2(h^2)}{dx^2} = 0$ (5.29)

如果应用裘布衣-福尔海默理论，把一个三维非封闭水流场降低为二维xy水平水流场，那么，均质、各向同性介质中水流的方程就变为

$\frac{\partial^2(h^2)}{\partial x^2} + \frac{\partial^2(h^2)}{\partial y^2} = 0$ (5.30)

换句话说，h²比h更能满足拉普拉斯方程。把稳定流边界值问题放在方裎(5.30)的基础上，用类比或数值模拟的方法去解决浅层、水平流场中的 h(x, y)值，就是可能的了。同时，对于非封闭含水层中的裘布衣自由表面水流来说，借把方程(2.77)左侧中的h代换为h²的方法，建立起一个非稳定水流方程，也是可能的了。

Harr(1962)详细研究了裘布衣-福尔海默(Dupuit-Forchheimer)理论的实践情况。Bear(1972)著作中包括有细致的理论处理。Kirkham(1967)考察了理论中的自相矛盾的问题并提出了一些揭示性的解释。这种方法被广泛地应用于工程实践中。

习题

考虑一个饱和、均质、各向词性、长方形、垂向剖面，ABCDA。其上部边界为AB，底部边界为DC，左侧边界为AD，右侧边界为BC，使DC的距离为AD的两倍。试对下列每一种情况绘出精确的定量流网。
1. BC和DC是不透水的。是一个不变水头边界，h = 100米。AD被分为两个相等的长度，其上半部分不透水，下半部分是不变水头边界，h=40米。
2. AD和BC是不透水的，AB是一个不变水头边界，h = 100米。DC被分为三个相等的长度，左侧和右侧部分不透水，中间部分为不变水头边界，h=40米。

使问题1中的垂向剖面ABCDA变为一个不规则四边形，方法是使B点沿垂向升起，而D和C点的高程为0米。A是100米，B是130米。使AD，DC，和BC是不透水的，并使AB代表一个不变坡度的地下水位(其上的水力水头等于其高程)。
1. 对于这种情况绘制出精确的定量流网并用它们正确的h值标注等势线。
2. 如果此区域内的水力传导率是0^-4米/秒，试计算通过系统的总水流(单位以立方米/秒表示)(垂直这个剖面单位厚度的)。
3. 在流线与上部边界交叉处，用达西定律计算每一交叉点处的入流或出流流速。

1. 对于均质、各向异性的情况重复问题2(a)和2(b)。其中水平水力传导率是0^-4米/秒。垂向水力传导率是0^-5米/秒。
2. 对(a)中的均质、各向异性层组画出水力传导率椭圆图。以椭圆上的适当结构表示出在你的流网中所指出的水流方向与水力梯度之间的关系，在流网上的两点是正确的。

对于如下情况重复问题2(a)，即自由表面水流(在大气压力下)排水点定位在BC的中点。排水是安排在与流网平面呈直角的方向上。

1. 对于双层的情况重复问题1(a)和1(b)，2(a)。其中流场的下面一半的水力传导率比上面一半的大5倍。
2. 对于双层的情况重复问题1(b),其中流场上面一半的水力传导率值比下面一半的大5倍。

在2、3、4、5题中的流网的流场中心处绘制一个压力计，并根据流网的情况给出压力计的水位。

1. 重新绘制图5.3的流网图，该图中大坝宽150米，位于一个厚度为120米的地层之上。h₁=150米，h₂=125米。
2. 如果地层为两层，上方60米的地层的渗透率为下方60米地层渗透率的十分之一。按7(a)重新绘制流网。

两个相距500米的压力计，底部分别位于一个非承压含水层的深度100米和120米的位置。在水平非渗透层之上，较浅的压力计的水位为170米，较深的压力计水位为150米。使用Dupuit-Forchheimer 假设来计算两个压力计中点的水位高度，并计算水流通过一个厚度为10米的区域的水量。该区域的水力传导系数为 K = 10^-3 m/s 。

绘制一个水平承压含水层的水平流网图：
1. 水流流向一个处于稳态的抽水井（比如，水位保持恒定的水井）。
2. 两个处于稳态，并且抽水流量相同的抽水井（比如，井中的水头相同）。
3. 水井是一个线性，并且水头恒定的边界。