图4.12 三角形单元的节点控制区
如图4.12所示，根据前述模型的建立过程，三角形单元所覆盖的区域分别属于三个顶点的水均衡控制区。在有限差分模型中，节点控制区由三角形边线的中垂线D—E—F围成；在有限元模型中，节点控制区由三角形边线中线与重心的连线A—B—C围成。
在有限差分模型中，流向节点i控制区的流量由Qj，i和Qk，i两部分组成，Qj，i是经过断面D—E流入的，而Qk，i是经过你断面E—F流入的。两部分流入的总和为
Qi?Qj,i?Qk,i?TDEEF(Hj?Hi)?T(Hk?Hi) ijik（4.3.62）
利用三角形中垂线的性质可以得到（陈崇希等，1989）
bbbb?ccDEij?cicjEF ??,??ikik （4.3.63） 4Ae4Aeijik式中：Ae为单元面积。将式（4.3.63）代入式（4.3.62）可得
Qi?T??bi(bj?bk)?ci(cj?ck)??4AeHi?T?bibj?cicj?4AeT?bibk?cick?Hj?Hk （4.3.64）
4Ae在有限元模型中，流向节点i控制区的流量只要把式（4.3.39）改写为以下形式即可
T(bbT(bi2?ci2)T(bbij?cicj)ik?cick)Qi??Hi?Hj?Hk
4Ae4Ae4Ae（4.3.65）
比较式（4.3.64）和式（4.3.65）可以发现，两者的后2项（关于Hj和Hk）完全相同
只有与Hi相乘的系数在形式上有所不同，但实际上也是相同的。利用式（4.3.12）和（4.3.13）给出的关系，很容易证明
bj?bk??bi,cj?ck??ci
（4.3.66）
因此，式（4.3.64）和式（4.3.65）给出的结果是完全一样的。这说明无论有限差分模型还是有限元模型，计算得到的侧流量相同，在这一点上，他们是等价的。然而他们又不是完全等价的，因此他们计算的控制区面积可能不同，只有在图4.12中点B和点E重合时，两者计
算的控制区面积才完全相同。
上述等价性质导致多边形网格有限差分模型的差分方程和三角形网格有限元模型的有限元方程具有相同的形式，都可以写为
kGppH??Gp,L?iHL?i?Qa?kpi?1NpSAp?tk?Hkpk?1?Hp?
（4.3.67）
式中：Gpp和Gpl是由式（4.3.44）得到的关联系数；L-i表示与节点p相邻的第i个节点；
Np是与p相邻的节点总数；Qa为其他外部补给量；Ap是节点水均衡控制区的面积；k为时步。需要注意的是，式（4.3.67）运用了储量集中处理的有限元法。对于关联系数，可以证明
kGpp???Gp,L?iHL?i
i?1Np（4.3.68）
当p点存在抽水井时，常规的有限元法是直接把井流加入到方程的Qa之中，同样导致“算不准”问题，需要校正。不妨把方程（4.3.67）中的Qa表示为
Qa??Qw??pAp （4.3.69） 式中：Qw为井孔排泄地下水的流量；?为其他面源的补给强度。进一步把方程（4.3.67）写成以下形式
其中：?因子定义为
kGppH??Gp,L?iHL?i?Qw??Ap
kpi?1Np（4.3.70）
???p?Skk?1(Hp?Hp) （4.3.71） ?tk引入井孔附近的解析解，即（4.2.54），可得
kkHL?H?iw?Qwr?2lnL?i?rL?i 2?Trw4T（4.3.72）
式中：rL-i为节点L-i与节点p的距离。将（4.3.72）代入式（4.3.69）得
Qw??ApB?A(H?H)?Qw?E?
TTTkpkw（4.3.73）
其中，A、B、E这三个系数分别为
A?GppT,B??i?1NpGp,L?irL?i1NPGp,L?i2ln,E??rL?i 2?Trw4i?1T（4.3.74）
可以证明，在井孔附近均质各向同性条件下，他们实际上与倒水系数T无关。这是因
为，关联系数在计算时已经包括了导水系数，见式（4.3.65）或式（4.3.40），当关联系数再
除以T时，倒水系数在A、B、E中被消去了。因此，这三个系数只是取决于模型网格几何信息的常数。如果引入图4.4中的wi和li表示网格的几何信息，Li=rL-i为节点L-i与节点p的距离，wi为两者连线的中垂线宽度，则有
于是上述三个集合因子可另外表示为
Gp,L?i?Twi rL?iN（4.3.75）
w1A??i,B?2?i?1rL?iNPwirL?i1pln(),E??(wirL?i) （4.3.76） ?rrw4i?1i?1L?iNp(a) n=3 （b）n=4 a (c)n=5 (d)n=6
图4.13 具有正n边形特征的节点控制区
由式（4.3.73）可以得到井口的校正公式
Qw?E?ApAkkT(Hp?Hw)?? B?1B?1（4.3.77）
当节点的水均衡控制区为正多边形时，几何因子E与控制区的面积Ap往往相等，从而在式
（4.3.77）中消去e因子的影响。图4.13给出了有限差分模型中节点控制区为正多边形的例子，设正多边形的定点数位n边长为a，则
wi?a,rL?i?atan()n?
（4.3.78）
?A?a12?A?ntan(),B?ln?,E?A?natan(?/n) （4.3.79） ?pn2??rwtan(?/n)?4校正公式简化为
其中：等效半径
Qw?kk2?T(Hp?Hw)ln(req/rw) （4.3.80）
req???2?exp?1??
tan(?/n)ntan(?/n)??a（4.3.81）
当n=3，4，5，6时分别得到
n?3,req?0.172a;n?4,req?0.208an?5,req?0.244a;n?6,req?0.282a
其中，当n=5和n=6时的等效半径适用于有限元模型。这些结果为在多边形网格有限差分
模型和有限元模型中模拟含水层—井孔系统水流问题提供了简单的校正方法。
五、控制条件刻画及方程组的求解
有限元模型对边界条件的处理基本上与基于多边行网格的有限差分法相同，在此不赘述。对于越流条件，把越流补给项加入到有限元方程中，则式（4.3.67）可更新为
kGppH??Gp,L?iHL?i?kpi?1NpKcApmc?Hc?Hpk??Qa?SAp?tk?Hkpk?1?Hp? （4.3.82）
式中：KC和mc分别为弱透水层的垂向渗透系数和厚度；Hc是作为补给源的含水层的水头；Ap为节点p控制区的面积。
把考虑各种条件的节点有限元方程汇总之后，可以形成形如式（4.3.67）的线性方程组，并可采用直接法或迭代法进行求解。
第四节 地下水流边界元模拟
边界元法（boundary element method）又称为边界积分方程法，字20世纪70年代提出之后，在地下水流数值模拟领域有所应用。其基本思路是利用格林（Green）公式和偏微分方程的格林函数，把地下水流问题转化为边界积分方程组的求解问题。 边界元法与有限差分法和有限元法的显著差异，是不必在模拟区域内部布置离散点，而只需要将模型的边界离散化，在求解出边界节点的水头和水力梯度之后，就可以用边界积分方程计算出区域内任一点的水头。
一 、格林公式与边界积分方程
对于在平面区域Ω内定义的一阶连续函数P（x，y）和Q（x，y），满足格林定律
??p?Q??Pcos?n,x??Qcos?n,x??ds ??dxdy?????????B??x?y? （4.4.1）
式中：B代表区域Ω的边界；n为边界的外法线方向。格林函数说明了平面二重积分与边界积分之间的关系。如果函数P（x，y）、Q（x，y）与一阶连续函数G（x，y）、H（x，y）之间存在如下关系
P?G?H?H;Q?G ?x?y（4.4.2）
则可以通过格林定律推导得到第一格林公式
??2
?G?2Hdxdy???GB??G?H?G?H??Hds?????dxdy ???n??x?x?y?y?（4.4.3）
式中：?为Laplace算子。另外还可以推导出第二格林公式
?G???H22??G?H?H?Gdxdy?G?Hds ??????B????n?n??（4.4.4）
这两个格林公式是建立边界元方法的基础。
设区域Ω内有均质各向同性的承压含水层，其地下水的稳定流偏微分方程为
??2H?2H?T?2?2????x,y??0
?y???x（4.4.5）
式中：??x,y?为源汇项。式（4.4.5）可进一步表示为
2
把式（4.4.6）代入第二格林公式有
?H????x,y?T （4.4.6）
???G??(x,y)T?G???Hdxdy???H?2Gdxdy??G?Hds ?B????n?n??（4.4.7）
进一步转化为
?H??(x,y)??G2 H?Gdxdy?H?Gds?Gdxdy （4.4.8）?????B??????n?nT??现在研究区域内固定点P的情况（图4.14），以该点为中心建立极坐标系，定义格林函数 它满足
G?x,y??G?r??lnr
?2G???r?
（4.4.9）
（4.4.10）
式中：??r?为Dirack函数。同时引入过渡函数
f?r,???H?2G?Hp??r?
（4.4.11）
式中：?为极坐标下的辐角;Hp为p点的水头于是可以得到
???H?Gdxdy??[lim?f(r,?)dr]d? （4.4.12）
0r?0022?r利用Dirack函数的性质，上述积分可计算为
2H????Gdxdy??2?0Hpd??2?Hp （4.4.13）