?1??xL，?2??xL，?3??xL，?4??xL （4.2.60）
1234其中：系数A、B、E都只取决于网格的形状，可称为几何因子。式（4.2.28）建立了井孔流量Qw、计算井点水头H0和井孔实际水头Hw之间的定量关系，即
Qw?AE??x?ykT(H0k?Hw)?? （4.2.61） B?1B?1这是综合考虑面源补给、地下水储存量变化和非正方形网格的条件下推导出来的校正公式，
因此更具有普遍意义。
下面讨论?因子的影响。如果井点单元的相邻单元都具有相同的横向尺寸寸?y，则有
和纵向尺
E??x?y （4.2.62）
因此，式（4.2.61）中关于?因子的项被抵消了，校正公式可进一步简化为
Qw?AkT(H0k?Hw) （4.2.63） B?1只要在模型网格剖分时稍加注意，这一点很容易实现。实际上，在相邻单元尺寸与井点单元尺寸相差不显著的情况下，几何因子E接近于井点单元的面积?x?y，从而使?因子对井孔校正结果的影响比较小，基本可以忽略。由于具有这种能够被抵消的性质，即使?因子是非均匀分布的，其影响也将在一定程度上削弱。当然，如果?因子的分布在井孔附近剧烈变化，假设?因子与径距无关还是可能带来一定的误差，有待改进。
当井孔附近的网格尺寸完全均匀时，?x??y，有
A?4，B?则校正公式可还原为
2?ln?x，E??x2 （4.2.64） rwk2?T(H0k?Hw)Qw? （4.2.65） reqln()rw其中等效半径req与式（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。这样就解决了常规有限差分模型的“算不准”问题。
关于多边形网格有限差分模型中井孔的校正，将在有限元模型的井孔校正方法中给出。
六、差分方程组的求解方法
设有限差分模型中又M个节点，对于任一节点p都可以建立形如式（4.2.26）或式（4.2.35）和式（4.2.36）的差分方程，不妨用下述形式的通式描述
kkkkk?1 （4.2.66） Cp0Hp?Cp1Hp?1?????CpiHp?i?????CpiHp?Np?Fp，p?1,2,???,M式中：下标p-i（注意中间并非减号，为两者连接号）为与节点p相邻的第i个节点在整体网格中的序号；Cpi为节点p和节点p-i之间的水流关联系数；Fp为对应节点p、取决于前一时刻水头分布和本时步之间步长的已知变量。考虑Cpi与待求变量Hp?i无关，则上述方程组为线性方程组，可用矩阵和向量表示为
??C???H???F? （4.2.67）
式中：??C??为渗透矩阵；为
?H?为节点水头列向量；?F?为右端项列向量。改方程组的解
?H????C???1?F? （4.2.68）
运用线性方程组的求解方法可得到上述解。
对于节点数M较小的问题，线性方程组的求解可采用Gauss消去法，LU分解法等直接法。然而，当节点数很大时，直接法的工作量太大，可考虑采用其他方法，如迭代法、共轭梯度法等。为了说明迭代法的思路，我们把式（4.2.66）中表示时刻的上标k暂时省略，方程组变为
Cp0Hp?Cp1Hp?1?????CpiHp?i?????CpiHp?Np?Fp，p?1,2,???,M （4.2.69）
关于节点p的方程可转化为
Hp?FpCp01p??CpiHp?i （4.2.70） Cp0i?1N上式提供了一种迭代求解的途径。迭代法的步骤是首先令节点的水头等于前一时刻的水头，即
Hp?Hp，p?1,2,???,M （4.2.71）
式中：上标0为第0次迭代运算。再调用式（4.2.70）进行迭代运算，每次迭代运算的公式为
Hn?1p0k?1?FpCp01pn?CpiHp??i，n?0,1,2??? （4.2.72） Cp0i?1N式中：上标n为第n次迭代。经过足够多次的迭代运算之后，前后两次迭代得到的水头值几乎相等，判断指标为
nn?1?Hmax?maxHp?Hp；p?1,2,???,M （4.2.73）
??式中：?Hmax为每次迭代结果与前一次迭代结果的最大差距，如果他小于允许误差，迭代
过程即可停止
?Hmax??H （4.2.74）
其中：?H是预先设置的收敛标准。
为了使迭代收敛更快，可以引入松弛因子?，并把迭代公式（4.2.72）改为
n?1Hp??(FpCp01pnn?CpiHp??i)?(1??)Hp，n?0,1,2??? （4.2.75） Cp0i?1N松弛因子的取值范围是1～2，这种迭代格式称为超松弛法。在采用迭代法求解某些非线性问题时，可能使用??1才能得到收敛结果，这种迭代格式称为低松弛法。
第三节 地下水流有限元模型
有限元数值模拟技术在固体力学和流体力学中都有广泛的应用，也经常用于进行地下水流的数值模拟。有限元法的原理与有限差分法有所不同，有限差分法以偏微分方程的差分近似为基础，而有限元法以积分方程的离散近似为基础。在历史上，曾经有两种推导有限元方程的途径：迦辽金（Galerkin）法与里茨（Ritz）法。Galerkin法采用加权余量的积分公式（以平面模型为例）为：
???'?R?H?(x，y)?w(x，y)dxdy?0 （4.3.1）
式中：R（H’）为近似解H’形成的余量函数；w为权函数；Ω为模型的某个子域。
Rietze法是利用与地下水偏微分方程等价的函数极限值函数进行求解，即
?H?HJ(H)???F(x，y，H，，)dxdy
??x?y （4.3.2）
式中：J(H)是等价范涵，积分内公式F定义为使J(H)取极小值的函数H（x，y，t）恰好满足地下水流偏微分方程。于是，地下水方程的近似解通过以下极值条件方程
?J?F???dxdy?0，p?1，2，3，??? （4.3.3）
??H?Hpp得到，其中Hp是模型子域上的节点水头。已有证明，当有限元网格单元采用相同的插值函
数及形函数时，上述两种方法实际上是完全等价的，所形成的代数方程组相同。为此，本书中不再对有限元法的积分方程进行推导。
一、单元与形函数
最常见的平面有限单元网格是三角形网格，而四边形网格也有较多的应用。三位有限元 模型一般采用分层网格建立6节点或8节点等参单元。有限元的积分方程就是在网格单元的基础上近似求解的。
yjkeiox
图4.9三角形单元及其顶点编号
首先讨论平面三角形单元（图4.9）。单元e的三个顶点编号按照逆时针顺序分别为i、j、k，水头在单元内任意坐标点（x，y）的分布特征采用线性插值函数进行近似描述
H(x，y)??1??2x??3y （4.3.4）
式中：β1、β2、β3为几何系数。把节点i、j、k处的坐标和水头代人（4.3.3）得
?1??2xi??3yi?Hi （4.3.5）
?1??2xj??3yj?Hj （4.3.6）
这是线性方程组，求解得
?1??2xk??3yk?Hk （4.3.7）
1(aiHi?ajHj?akHk) （4.3.8） 2Ae1(biHi?bjHj?bkHk) （4.3.9） 2Ae1(ciHi?cjHj?ckHk) （4.3.10） 2Ae?1??2??3?其中：
ai?xjyk?xkyj；aj?xkyi?xiyk；ak?xiyj?xjyi （4.3.11）
bi?yj?yk；bj?yk?yi；bk?yi?yj （4.3.12） ci?xk?xj；cj?xi?xk；ck?xj?xi （4.3.13）
Ae?(bcij?cibj) （4.3.14）
式中：Ae为三角形单元面积。将式（4.3.8）至式（4.3.10）代人式（4.3.4），则插值函数可
改写为
H(x，y)?HiNi(x，y)?HjNj(x，y)?HkNk(x，y) （4.3.15）
其中:
Np(x，y)?1(ap?bpx?cpy)，p?i，j，k （4.3.16） 2Ae是节点i、j、k在单元e内的形参数。对于固定坐标点（x，y）形函数具有以下的性质
eNie?Ne?Njk?1 （4.3.17）
式中：上标e表示求和只能在同一三角形单元内进行。当所有节点的水头已知时，就可以把
形参函数代人式（4.3.15）计算三角形单元内任意坐标位置的水头。
y12ξξ=141η=-1o4ξ=-1η=123ηo(a)x3(b)
图4.10任意四边形单元的坐标变换
四边形单元也可以采用类似的方法建立插值函数和形函数，但是在形式上复杂一些。首先，任意的四边形都可以影射为一个正方形，图4.10（a）中（x，y）坐标系下的四边形1-2-3-4，可通过变换，形成4.10（b）中(?，?)坐标系下的正方形1-2-3-4。这个正方形单元被称为等参单元，其形函数定义为
Ni(?，?)?1(1???，2，3，4 （4.3.18） i)(1??i?)，i?14式中：(?i，?i)为节点i的影射坐标。等参单元内的任意一点(?，?)，通过上述形参数与实际四边形单元内的对应坐标点建立如下关系，即
x??xiNi(?，?)，y??yiNi(?，?) （4.3.19）
i?1i?144式中：xi、yi为节点i的实际坐标。因此，等参单元内形函数与实际单元函数具有一一对应的关系
Ni(?，?)?Ni(x，y) （4.3.20）