孙婧祎,周德亮,张思宇
(辽宁师范大学 数学学院,辽宁 大连 116029)
井水位计算在地下水工作中具有重要的作用,目前已有多种计算井水位的方法。一些从事地下水数值模拟的人员直接将传统方法求得的节点水位看做是该点的井水位。这是不正确的[1]。因为传统方法计算得到的井水位只能代表井点所在节点均衡区域上的平均水位[2],并非井的真实水位,所以应修正传统地下水模拟中对井的处理。在已经产生的多种无网格方法中,径向基函数配点法[3][4]具有易实施、高精度、高效率、与空间维数无关、不受背景网格限制等优点,是一种真正的无网格方法。本文的工作是将径向基函数配点法用于求解承压稳定流井水位计算的问题中,并对传统方法、传统方法的校正方法[6]及奇点磨光法[7]这三种计算井水位的方法进行了比较。
具体考虑如下平面承压稳定井流混合问题[5]
(1)
该问题中依次出现的函数及变量的含义为:T是导水系数,h(x,y)是给定的水头函数,ε(x,y)和q分别表示为越流补给强度和补给量,qk和ρk则分别表示第口井的开采量和井口半径的大小。
(2)
在上式中Cj(j=1,2,…,N)是待求系数,令x=(x,y),xj=(xj,yj)。
(3)
(4)
(5)
进而,可写成线性方程组的形式,即
Ac=b
(6)
其中A=[Lφ L2φ φ],b为(5)式中方程的右端项。
由(6)式得到方程组的系数,代入到(4)式中即可得到该混合问题的数值解。
传统方法是将问题(1)转化为如下的等价问题,即
(7)
实际水流是径向入流,水头呈非线性分布。因此,用传统方法模拟得到的井点水位h,并不能代表井中的真实水位hw,而是相当于距离井点rf为处的实际水头,如图1[6]。这里需要注意的是,本文所讨论的仅是矩形网格中井水位校正的问题。
图1 井点附近水头剖面分布
对于井点附近划分的网格间距恒为Δx,Δy的矩形网格,根据Theim公式,井中水位hw与节点j,k,l,m(包围井点所在节点的四周节点)分别对应的水位hj,hk,hl,hm之间满足如下关系,即
(8)
(9)
其中,Q为抽水井的流量,rw为抽水井的有效半径,T为含水层导水系数。
数值法模拟得到i点水位为h,并且从周围单元进入i点单元的流量与井流量相同,那么由达西定律可得
(10)
把(8、9)式代入到(10)式中,求出hw,得到如下校正公式
(11)
(12)
由此便求得w,进而得到问题(1)的解h,即h=w+v。
图2 计算域及节点配置图
图3 井水位计算的传统方法的近似解 图4 传统方法校正方法的近似解
图5 奇点磨光法的近似解 图6 解析解
(13)
本文所列举的三种方法虽然可以满足一般地下水模拟与预测中对井水位计算的需要,但传统方法在井点处的计算误差较大,并且为了计算井点的控制域面积,就需在井点附近进行局部的网格剖分,这破坏了无网格方法的特性;传统方法的校正方法在计算井水位时可以达到很好的模拟效果,但在求解精度上低于奇点磨光法,另外需要注意的是本文所说的校正方法还仅是针对在井点附近进行矩形网格剖分的情况。因此,结合以上叙述并针对本文所做的实例,在基于径向基函数配点法计算井水位的问题中,采用奇点磨光法效果最好。该数值例子也说明了奇点磨光法的有效性。本文所给出的三种井水位计算方法也可应用到径向基函数配点法中的其它问题,这将在以后的工作中进一步讨论。
[1]陈崇希,唐仲华.地下水流动问题数值方法[M].武汉:中国地质大学出版社.1990:239.
[2]薛禹群,谢春红.地下水数值模拟[M].北京:科学出版社.2007:246.
[3]Sharan M, Kansa E J, Gupta S. Application of multiquadric method for numerical solution of elliptic partial differential equations[J]. Appl Math Comput, 1997, 84: 275-303.
[4]Franke C, Schaback R. Solving partial differential equations by collocation using radial basis functions[J]. Appl Math Comput, 1998, 93: 73-82.
[5]孙讷正.地下水流的数学模型和数值方法[M]. 北京:地质出版社.1981.7-48,110.
[6]陈崇希,王旭升,胡立堂.地下水数值模拟中抽水井水位校正[J].水利学报.2007,38(4):481-485.
[7]吉林大学数学系计算服务站地下水小组.地下水非稳定流计算的有限元方法(I)[J].吉林大学报;自然科学版.1977(1):15-27.