而等参单元内的任意一点(?，?)的水头，也可以通过形函数 变换为坐标下(x，y)处的水头
H(?，?)?H(x，y)??HiNi(?，?)??HiNi(x，y) （4.3.21）
i?1i?144在建立有限元方程时，需要用到一下类型的偏导数
4?N(???)?f （4.3.22） ??fii?xi?1?x?Ni(???)?Ni???Ni???? （4.3.23） ?x???x???x?Ni1?N1??i(1??i?)，i??i(1???（4.3.24） i)??4??4
二、基于三角形单元的地下水流有限方程
对于三角形网格中的一个节点p（图4.11），建立有限元方程的积分区域包括了与p点相邻的所有节点p-1、p-2等，即以p点为共同顶点的所有三角形单元所覆盖的区域。积分采用分片发，即在相邻的三角形单元内分别积分，然后求总和，得到对应于p点的有限元方程。
p-4e-3yp-3pe-4e-2p-5p-2CABVe-1x p-2e-5p-1e-1p-1o图4.11三角形网格中的节点及其水均衡控制区（阴影部分）
通过建立有限元网格节点的水均衡控制区，采用水量均衡原理，也可推导出地下水流的有限元方程，与积分法得到的结果相同。这种推导方法更加清楚的反映了地下水流的物理机制。如图4.11所示，节点p的谁均衡控制区是一系列三角形重心与侧边中点的连线所围成的区域，每个三角形单元中，属于节点p谁均衡控制区的面积是其单元面积的1/3。因此，节点p谁均衡控制区的面积为
1p（4.3.25） Ap??Ae?n3n?1
式中：e-n为第n个相邻三角形单元的整体编号；Ae-n为其面积。 地下水通过控制区周边流入控制区内侧向流量QH，可以用流速在控制区边界线上的线积分得到
pQH?m??S(V?n)dS?m??S??Vxcos(n，x)?Vycos(n，y)??dS （4.3.26）
N式中：V为流速向量；n为边界线S上微分线段dS的内法线矢量；Vx和Vy分别为x方向和y方向的流速分量；m为含水层的厚度。设渗透系数的主轴与坐标轴一致，则流速向量可根据达西定律获取
Vx??Kx代人式（4.3.26）有
pQH??m??[KxS?H?H，Vy??Ky（4.3.27） ?x?y
?H?Hcos(n,x)]dS?m?[Kcos(n,y)]dS（4.3.28） y?S?x?y
Np确定上述积分的方法是在节点p相邻的三角形单元内分别计算，然后求和，即
e?nQ??QHpHn?1（4.3.29）
其中，每个三角形单元e-n流量贡献为
e?nQH??[Kxm???H?H|e?ncos(n,x)dS]?[Kym??S?n?y|e?ncos(n,y)dS]（4.3.30） S?n?x式中：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），有
?N?N?N?H?Hpp?Hp?1p?1?Hp?2p?2（4.3.31） ?x?x?x?x
根据式（4.3.16）得
?H1（4.3.32） ?(Hpbp?Hp?1bp?1?Hp?2bp?2)?x2Ae?1
同理有
?H1（4.3.33） ?(Hpcp?Hp?1cp?1?Hp?2cp?2)?y2Ae?1
可见水力梯度在三角形单元内是一常量，可以提取到式（4.3.30）的积分号外部，即
?Km??Kym?e?1QH???x(Hpbp?Hp?1bp?1?Hp?2bp?2)?cos(n，x)dS?(Hc?Hc?Hc)cos(n，y)dS???ppp?1p?1p?2p?2???S?1S?12A2A?e?1??e?1?（4.3.34）
这样在线积分中只剩下积分路径在x轴和y轴上的投影，有
??S?1cos(n，x)dS?yA?yC，??cos(n，y)dS?xC?xA （4.3.35）
S?1根据点A和点Ｃ与单元e-1的三个顶点之间的关系，有
11(xp?xp?1)，yA?(yp?yp?1) （4.3.36） 2211xC?(xp?xp?2)，yC?(yp?yp?2) （4.3.37）
22xA?将式（4.3.36）和式（4.3.37）代人式（4.3.35），得
??2Txbp?Tymc2pS?1cos(n，x)dS?cp，??cos(n，y)dS??bp （4.3.38）
S?1
将式（4.3.38）代人时（4.3.34），得
Qe?1H=
?4Ae?1Hp?Txbpbp?1?Tycpcp?14Ae?1Hp?1?Txbpbp?2?Tycpcp?24Ae?1Hp?2（4.3.39）
三角形单元e-1中的节点p、节点p-1以及节点p-2分别与图4.9中普遍单元的顶点i、j、k相互对应，有了式（4.3.39）之后，引入如下的系数来反映单元e内节点与节点之间的关系
eGPL??eeeTxbpbL?TycepcL4Ae，p?i,j,k，L?i,j,k （4.3.40）
式中：Cpl为节点L与节点p的控制区之间得到关系。利用式（4.3.40）可以把式（4.3.39）
改写为
e?1e?1e?1e?1QH?Gp,pHp?Gp,p?1Hp?1?Gp,p?2Hp?2 （4.3.41）
代人式（4.3.39）得到流向节点p控制区的总侧向流量
pe?ne?ne?nGH??(Gp,pHp?Gp,p?jHp?j?Gp,p?kHp?k)n?1Np（4.3.42）
式中：p-j和pk分别表示沿着单元e-i从节点p逆时针移动遇到的第一个节点与第二个节点。
以图4.11中的节点p为例，展开式（4.3.42），得到
pe?ne?1e?5e?1e?2e?2e?3GH?(?Gp,p)Hp?(Gp,p?1?Gp,p?1)Hp?1?(Gp,p?2?Gp,p?2)Hp?2?(Gp,p?3?Gp,p?3)Hp?3n?1e?3e?4e?4e?5?(Gp,p?4?Gp,p?4)Hp?4?(Gp,p?5?Gp,p?5)Hp?5Np（4.3.43）
在整个网格中，节点之间在各个单元中的关联系数可以合并为
eGpL??GpLe?1EpL （4.3.44）
式中：Gpl为节点p与节点L之间在整体网络上的关联系数；Epl为共同拥有节点p和L的单元的数目。，说明节点p和L之间没有直接关系。关联系数的重要特征是对称性：Epl=Gpl。运用式（4.3.44），流向节点p控制区的侧向流量可表示为
Q??GpLHLpHL?1N（4.3.45）
式中：N为网格中的总节点数。 对于承压含水层，在时间步长#tk内，节点p控制区内水体积的增加量为
kk?1??Vw???S?H?Hx,yx,y?dxdy （4.3.46） ??式中：k为时间步长。标准的有限元法是把差值函数引入上述积分，并分别在相邻单元内积
分再求和，导出#Vw中包含全部相邻节点的水头变化信息。目前，地下水流有限元模型中，普遍采用储量集中处理法，即用节点p的水头变化表示控制区的平均水头变化，有
kk?1?Vw?SAp(Hp?Hp) （4.3.47）
其中：Ap按照式（4.3.25）计算，这种简单的做法反而可以改善有限元的计算结果。根据书均衡原理，控制区水体积的增量与外部流量补给的水量相等，即
P(QH?Qa)?tk??Vw （4.3.48）
式中：Qa是其他外部不计流量。利用式（4.3.45）和式（4.3.47），可以展开式（4.3.48）为
(SAp?tk?Gpp)H?kpL?1,L?p?NkGpLHL?Qa?SAp?tkk?1Hp（4.3.49）
这就是承压含水层模型中对应节点p的隐式有限元方程。 潜水含水层的饮食有限元方程课写为
(SyAp?tk?Gpp)H?kpL?1,L?p?NGpLH?Qa?kLSyAp?tkk?1Hp（4.3.50）
式中：Sy为给水度。关联系数GpL中的导水系数Tx、Ty与当前时刻的带球水头有关，因此属于非线性方程。
三、三角形四边形混合单元地下水流有限元方程
对于三角形单元，已经通过式（4.3.40）得到了单元内e内关联系数Gpl的表达式。在有限元模型中，对于任意单元体，节点之间关联系数的积分形式为
eGpLeeee???Np?Np?NL?NL????Tx?Tydxdy （4.3.51） ????x?x?y?y?e?
这个积分表达式对三角形单元同样是有效的。对于四边形单元，需要将上述积分转换到射影
的等单元上进行，利用以下性质
dxdy?式中：
Jd?d? （4.3.52）
J为坐标变换的Jacob矩阵，而
J为其行列式，有
??x/???? ?J?=???x/??并且
?y/????
?y/??? （4.3.53）
ee??NL?1??NL/???/?x???J? ? （4.3.54）
??Ne/?y???????Ne/?????L??L??11??y/???? J= ???J???x/????y/???? （4.3.55）
?x/???而x、y对η、ε的偏导数具有式（4.3.22）的形式。
利用式（4.3.52），积分式（4.3.51）可转变为
GepLeeee???Np?Np?NL?NL????Tx?Ty????x?x?y?y?1?1??11Jd?d?
（4.3.56）
注意积分符号内的偏导数和
J也是影射坐标（η，ε）的函数。
上述积分可采用高斯积分公式进行计算
??11?1?1f(???)d?d??i?1,j?1?2,2f(?i??j)
（4.3.57）
其中ηi，εj），i=1，2，3，4为四个高斯积分点，并有
?1??1??11；?2??2? 33（4.3.58）
在确定了单元内的节点关联系数之后，就可以利用（4.3.44）组装整体关联系数。而
出水量部分仍然采用储量集中处理法，节点p控制区在单元e内所占的面积为
A???Ndxdy??eepep1?1?1?1eNp|J|d?d? （4.3.59）
而节点p控制区的总面积是
其中：Ne,p为与p相邻的单元数。
最后可以得到承压水层三角形四边形混合单元模型的有限元方程为
e Ap??Apn?1Ne,p（4.3.60）
?GpLHLk?Qa?l?1NSAp?tk?Hkpk?1?Hp?,p?1,2,...,N
（4.3.61）
四、地下水流有限元模型的孔井校正
与有限差分模型一样，地下水流有限元模型在遇到井孔时也需要进行校正。井孔的校正思路与有限差分模型中所采取的思路是相同的．即引入井孔附近的径向流解析解．在研究地下水流有限差分模型的井孔校正问题时，我们已经建立厂矩形网格有限差分模型的校正公式，但是尚未对多边形网格有限差分模型进行推导。在本节中．我们将证明，在各向同性条件下．多边形网格有限差分模型与三角形网格有限元模塑在某种程度上是等价的，因此，两者的井孔校正问题可以一起处理。