地下水运动的数学模型

第四章 地下水运动的数值模型

解析解虽然具有精确可靠的特点，但采用解析解反映自然状态和复杂人类活动干扰下的地下水运动是相当困难的。因此，当含水层的条件严重偏离现有解析模型的简化假设时，人们通过数值模型来获得近似的地下水流场及演变趋势。

第一节 地下水流数值方法概述

地下水流的数学模型采用偏微分方程描述地下水流的时间和空间连续状态，而数值模型则是采用离散（非连续）时空模型中水头的分布与演变对数学模型进行近似描述。从精确数学模型到近似数值模型的转化，虽然会损失一些精度，但使复杂地下水流问题的分析得以通过机械计算实现，而且误差也是可控的。

把偏微分方程求解的数值方法引入到地下水流问题的求解始于20世纪70年代，主要方法包括有限差分法、有限元法和边界元法，此后又发展了有限分析法、多重网格法和无网格法等。这些方法的共同特点是将模型空间及边界离散为由一系列的节点以及联系这些节点的单元（无网格法除外），含水层的水头在这些节点上定义，从而实现了水头分布空间连续函数向离散变量的转化，表示为

21

212111

22111221

202()

02()02()002(0)k k k k k k k k k k k k k k k k k k k k k k k k

k

k f f f f a b c e x L x x t t t t f x f f f f a b c e x L x x t t f f f f a b c x L e x x

d f df

e e

f a b f c x L dx dx t t f x u

---------∂∂-++=<<∂∂∆∆=-∂∂-++=<<∂∂∆∆∂∂=+

++<<∂∂+-++=<<∆∆==，，，，

{}(,,);1,2,3,,p H x y z H p M ⇒=⋅⋅⋅ （4.1.1）

式中；H 为含水层的水头；x 、y 、z 为空间坐标；p 为数值模型的节点；M 为节点的数目。与此同时，时间也被离散化为一系列具有先后顺序的时刻，不同时刻节点的水头分布方式也不同，每个时刻节点的水头分布结果总是由前一个时刻演变过来，即 {}{}

1

;1,2,3,,;;;1,2,3,,;;k k k k k p p p M t t p M t t t H H +=⋅⋅⋅=→=⋅⋅⋅=+∆ （4.1.2） 式中：k 为时刻；t 为时间；k t ∆为时间步长。数值模型的核心部分，是在地下水流数学模型的基础上确定节点在流场中的相互关系，并形成如下形式的方程组

11111212,1,2,k k k k k

p M p p p pp pM p M C C C C H H H H F +++++++⋅⋅⋅+==⋅⋅⋅⋅⋅⋅+（4.1.3） 式中：系数pL C 反映了节点P 与节点L 在流场中的关系，可称为关联系数；k p F 通常是由前一个时刻决定的常量。式（4.1.3）除包括地下水流控制方程的离散近似以外，还包括了边界条件的离散化处理。当k=0时，k p F 中包含了初始条件。如果模拟稳定流，则k p F 与时间无关，上标k 可以略去。使用矩阵符号，式（4.1.3）可缩写为

[]{}{}M M C H F ⨯= （4.1.4）

实际在处理模型时，一个节点通常只与相隔最近的若干节点有直接的关系，在式（4.1.3）中大多数关联系数为零，而且一般有对称性，即PL LP C C =。因此，矩阵[]C 通常是对称的

稀疏矩阵，降低了方程组求解的困难。尽管如此，区域地下水流数值模型的节点数往往很大，而且在某些非线性条件下，[]C 本身是待求水头的函数，方程组的求解需要使用迭代法而耗费相当多的时间。

不同类型的地下水流数值模型，最核心的差异在于形成式（4.1.3）的方式，关联系数具有不同的计算公式。原则上，对于完全相同的地下水流数学模型，以及完全相同的离散节点序列和离散时间序列，不同的数值模型应给出近乎同样的结果，只是计算精度略有差异。已有研究表明，某种类型的有限差分法与某种类型的有限元法存在一定的等价性，但更多的对比分析尚待进一步探讨。实际上，研究者在检验数值模型的可靠性时，总是与解析解的模型对比。当不同方法的数值模型用相同的解析解作验证时，他们之间的等价性就可以得到一定程度的证明。然而，人们对数值模型的偏好并不取决于数值模型之间的等价性，而往往取决于某种数值方法的熟悉程度，或者基于某种数值方法专业软件的方便快捷。

专业软件在地下水流数值模型的推广应用中发挥了重要作用，如有限差分模拟程序MODFLOW （McDonald and Harbaugh ，1988）和HST3D （Kipp ，1987）、有限元模拟程序FEMWATER （Lin et al ．，1997）和FEFLOW （Diersch and Kolditz ，1998）等，国内有陈崇希等开发的有限差分模拟程序PGMS （陈崇希等，2007）。其中，以美国地质调查局发布的有限差分模拟程序MODFLOW 最为著名，本章最后将介绍其特点。

第二节 地下水流有限差分模型

有限差分法是地下水流数值模拟最早使用的一种方法。这种方法古典而简洁、物理意义明了、编程容易，因而被广泛流传。地下水流的有限差分模型已经从最初的一维均质模型发展为目前的三维非均质模型，并有专业软件加以实现。

一、有限差分法的基本原理

有限差分法是求解偏微分方程边值问题和初值问题的一种数值方法，其实质是把连续的模型空间离散化为规则或不规则的网格点，利用导数的差分近似形式代替偏微分方程形成差分方程组，通过求解方程组得到离散点的待求变量作为连续场的一种近似结果。

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

首先用一维问题来阐述导数的差分格式。设待求变量f 沿x 轴的分布为二阶连续函数()f x ，现将x 轴离散化为节点012(,,,,,,)i N x x x x x ⋅⋅⋅⋅⋅⋅，它们的待求变量为012(,,,,,,)i N f f f f f ⋅⋅⋅⋅⋅⋅，则()f x 在节点i x 处的一阶导数表示为