把（4.4.13）带入（4.4.8）可得 2?HP=??[HB?G?H?(x,y)-G]ds???Gdxdy （4.4.14）
??n?nT这就是内点的积分方程。
ΒΒΩrα(x,y)pr
图4.14 模型区域内点和边界点
当p点位于区域边界上时（图4.14），如果边界是光滑的，则式（4.4.13）中的积分范围只能是边界切线的内侧，因此
pΩ(x,y)?HP=??[HB?G?H?(x,y)-G]ds???Gdxdy （4.4.15）
??n?nT?G?H?(x,y)-G]ds???Gdxdy （4.4.16）
??n?nT如果p点刚好是一个折线交点，内夹角为?，则
?pHP???B[H上式即为边界点的积分方程。
二、地下水稳定流的边界元方程
将模型区域离散化为一系列边界节点和线段，外边界顺时针编号，内边界逆时针编号（图4.15）。
isi+1内边界14131112p23
1图4.15 边界元模型的离散化方式
其中，每两个节点之间的线段即为边界元。边界节点的性质有两种：已知水头边界；已知水力梯度边界。令
J??H?G,U? （4.4.17） ?n?n则任意边界节点p的积分方程可离散化表示为
?pHP???e?1NB?e(HUP?GPJ)ds???Gp??(x,y)T dxdy （4.4.18）
式中：e为边界元的序号；N为边界元的书目；B-e为边界元e对应的线段。
设边界元e的起点为i，终点为i+1，它们的距离为Di,i?1。沿着边界元从i点向i+1点方向移动的距离定义为s，则可以利用向量几何关系建立以下函数（水头和水力梯度采用线
性插值）
H(s)?Hi?Hi?1?HiJ?Js, J(s)?Ji?i?1is （4.4.19）
Di,i?1Di,i?11?GLGP(s)?lnr2， UP(s)??2 （4.4.20）
2?nrr2?(s?s0)2?L2，Di,j?1?(xi?1?xi)2?(yi?1?yi)2 （4.4.21） L?xi?1?xiy?y(yi?yp)?i?1i(xi?xp) （4.4.22） Di,i?1Di,i?1xi?1?xiy?y(xi?xp)?i?1i(yi?yp) （4.4.23） Di,i?1Di,i?1s0?这些函数关系式用来计算（4.4.18）中的积分
?B-e(HUP?GPJ)ds??Di,i?1Di,i?10H(s)UP(s)ds??Di,i?10GP(s)J(s) （4.4.24）
把式（4.4.19）至式（4.4.23）带入式（4.4.24），得到积分的具体表达式为
?B-e(HUP?GPJ)ds??0(Hi?Di,i?11Hi?1?HiJi?1?JiL22s)ds?ln[(s?s)?L](J?s)ds0i22?0Di,i?1(s?s0)?L2Di,i?1 （4.4.25）
引入以下定积分记号
IB?e11LDi,i?1sds （4.4.26） ?Di,i?1?0(s?s0)2?L2LDi,i?1ds （4.4.27） 22?0Di,i?1(s?s0)?LB?eI12?IB?e21LDi,i?1s?ln[(s?s0)2?L2]ds （4.4.28） ?Di,i?102LDi,i?1122ln[(s?s)?L]ds （4.4.29） 0?0Di,i?12B?eI22?则式（4.4.25）可以转变为
?B-eB?eB?eB?eB?eB?eB?e(HUP?GPJ)ds?(Di,i?1I12?I11)Hi?I11Hi+1?(Di,i?1I22-I21)Ji?I21Ji?1
（4.4.30）
其中，4个定积分的公式如下
B?eI11(Di,i?1?s0)2?L2LB?e?ln,L?0;I11?0,L?0 （4.4.31） 222Di,i?1s0?LDi,i?1?s01sB?e?[arctan()?arctan(0)],L?0;I12?0,L?0 （4.4.32） Di,i?1LL
IB?e12B?eI21?[(Di,i?1?s0)2?L2]{ln[(Di,i?1?s0)2?L2]?1}?(s02?L2)[ln(s02?L2)?1]4Di,i?1(Di,i?1?s0)ln[(Di,i?1?s0)2?L2]?s0ln(s02?L2)2Di,i?1 （4.4.33）
IB?e22?B?e?I12L?1 （4.4.34）
当B-e中包含节点p时，注意运用以下极限 limxlnx?0
x?0把式（4.4.30）代入式（4.4.18）得
B?eB?eB?eB?eB?eB?e?PHP??[(Di,i?1I12?I11)Hi?I11Hi?1?(Di,i?1I22?I21)Ji?I21Ji?1]???Gpe?1?N?(x,y)Tdxdy （4.4.35） 这就是地下水稳定流对应节点p的边界元方程。对每个节点p都可以建立一个边界元方程，该节点的待求变量或为HP，或为JP，因此待求变量总数与边界元方程的总数相等，有唯一解。
三、地下水非稳定流的边界元方程
设区域?内有均质各向同性的第一类越流系统含水层，其他地下水的非稳定流偏微分方程为
?2H?2H?(x,y)Hc?HS?H （4.4.36） (2?2)???2?x?yTBT?t式中：?(x,y)为源汇项；B为越流因子。
对于这种非稳定流问题，可以用Laplace变换法解决。其思路是：将非稳定流偏微分方
程通过Laplace变换消去时间偏导数，然后采用稳定流边界元方程求解，但所适用的Green函数与式（4.4.9）不同，最后通过Laplace逆变换得到非稳定流的近似解。然而，Laplace变换的逆变换比较复杂，削弱了这种方法的实用性。
另外一种方法是先对时间偏导数进行差分近似，再采用类似稳定流的边界元模型分时步计算离散时刻的水头分布。采用隐式差分格式，方程（4.4.36）变为
?2Hk?2Hk?(x,y)Hc?HkSHk?Hk?1 （4.4.37） （2+）+??22?x?yTBT?tk在H中略去上标，并把上式进一步转化为
2?2H?ChH?k?(x,y)T （4.4.38）
其中：
THc1SHk?1 （4.4.39） Ch??,?(x,y)??(x,y)?2?S2BTB?tk把式（4.4.38）代入第二Green公式，得到
??f(r,?)
?2H(?2G?ChG)dxdy???[HB?H?H?(x,y)?G]ds???Gdxdy （4.4.40）
??n?nT对研究区域内固定点p，以该点为中心建立极坐标系，引入Green函数G及过渡函数
G(x,y)?G(r)??K0(Chr), 2f(r,?)?H(?2G?ChG) （4.4.41）
式中：K0为零阶第二类变形Bessel函数。可以证明过渡函数为
f(r,?)?Hp?(r) 于是有
（4.4.42）
???2H(?2G?ChG)dxdy??[lim?HP?(r)dr]d??2?HP （4.4.43）
0r?002?r将其代入（4.4.41）得
2?HP???[HB?H?H?(x,y)?G]ds???Gdxdy （4.4.44）
??n?nT这一方程与式（4.4.14）在形式上相同，但采用的Green函数G不同。同理可得形式上相同
的边界积分方程。
下面采用与图4.15相同的方式建立非稳定流的离散元模型。定义以下函数（水头和水力梯度采用线性插值）
H(s)?Hi?Hi?1?HiJ?Js, J(s)?Ji?i?1is （4.4.45）
Di,i?1Di,i?1?GL?ChK1(Chr) （4.4.46） ?nrGp(s)??K0(Chr), UP(s)?而r、s0、L的计算式与式（4.4.21）至式（4.4.23）相同。引入以下积分记号
IB?e11ChLDi,i?1sK1(Chr)ds （4.4.47） ?22Di,i?1?0(s?s0)?LChLDi,i?1K1(Chr)ds （4.4.48） ?220Di,i?1(s?s0)?L1Di,i?122sK[C(s?s)?L]ds （4.4.49） 0h0?0Di,i?1B?eI12?B?eI21?IB?e221Di,i?1?K0[Ch(s?s0)2?L2]ds （4.4.50） ?Di,i?10则节点p的边界元方程可写为
B?eB?eB?eB?eB?eB?e?pHP??[(Di,i?1I12?I11)Hi?I11Hi?1?(Di,i?1I22?I21)Ji?I21Ji?1]???Gpe?1?N?(x,y)Tdxdy（4.4.51）
需要说明的是，方程中Ji和Ji?1的系数项符号与式（4.4.35）相反。对每个节点p都可以建立一个边界元方程，该节点的待求变量或为HP，或为Jp，因此待求变量总数与边界元方程的总数相等，有唯一解。在转入下一个时步时，边界积分在方程中的各项系数并不变，但函数?(x,y)将发生变化。
四、源汇项的处理