孙玉平,周德亮
(辽宁师范大学数学学院,辽宁 大连116029)
目前计算地下水水头函数的主流方法有有限元和有限差分法,用它们解决问题时,都要把计算域离散为有限单元或网格,最后求出水头函数。径向基函数配点法是用径向基函数构造插值函数并采用配点法作为离散方案的一种无网格方法,与空间维数无关,此方法不需要背景网格,效率高,形式简单。本文利用的最小二乘配点法,是在径向基函数配点法的基础上,采用最小二乘法求解计算改进而来,对计算域进行节点离散,并布置辅助点,这样计算结果比径向基配点法精度更高。
在配点型无网格法中,微分方程只在计算域内各节点处严格满足,因此可能会产生较大误差。张雄[1]等人除了节点又在计算域内引入了辅助点。近似函数仍然只通过节点构造,但要求微分方程在所有节点和辅助点上满足。此时方程个数大于未知数个数,需要用最小二乘方法求解,因此将此方法称为最小二乘配点无网格法。
假设求解一个一维地下水稳定流问题,其对应渗流微分方程为:
式中:h为水头,K为渗透系数,ε为垂直向补给,Ω为研究区,g(x)为给定的函数,Γ为研究区Ω的边界。
对研究区 Ω 布置节点 x1,x2,…,xN,其中 x1,x2,…,xN0为内部节点,xN0+1,xN0+2,…,xN0+M,为 M 个辅助点,xN0+M+1,xN0+M+2为一类边界节点,总节点数为N=N0+M+2。
令所有节点满足微分方程(1),边界上的节点满足边界条件(2),
式(3)中有N0+M+2个方程,N0个未知数,需要用最小二乘法求解。为了保证解的精度,边界条件必须严格满足,因此将上式改写成
式(4)中h1为边界Γ上的节点参数组成的向量,h2为域Ω内的N0+M个节点的参数组成的向量。上式可改写为
式(5)对应于边界条件,(6)式对应于微分方程,令边界条件严格满足,由(5)可解得
将(7)式代入(6)式中,得
式(8)中有N0个未知数,N0+M个方程,因此没有通常意义下的解,需要在最小二乘意义下求解。在Matlab中用矩阵除法计算此类方程组,会自动利用最小二乘法求解,因此可以求出水头h(x)的近似解。而径向基函数配点法,只是在研究域内不添加辅助点。
自然界的多孔介质很多都是非承压的,最小二乘配点法可以很好地用来解决非承压含水层中的水流问题。下面通过此方法在两种非承压条件下的应用来说明这种方法的优越性。
如下图1[2]所示,假设有一非承压含水层,底板是水平不透水层,有降雨补给ε=0.004 m/d,含水层两端有两条河流流过,两河相距5 000 m,上、下边界分别为断面线Ⅰ、Π,水头分别为 h1、h2,设 h1> h2,h1=280 m、h2=220 m,土壤渗透系数K=0.785 m/d,在计算区域内满足
分别使用解析方法、径向基函数配点法、最小二乘配点法求解。图2为不同方法的计算结果,从图中可以看出,最小二乘配点法比径向基函数配点法有更好的精度,更为接近解析解的结果。
图1
图2
如图3[2]所示,设有一水库,初始状态下库岸地下水位和库水位齐平,使库水位由骤降至,这样库岸内地下水向库内排泄,潜水面形成一向水库降落的曲线,随着时间的增加,曲线向岸坡远处延伸,直至最终稳定下来。在区域内满足
其中土壤的渗透系数K=0.864 m/d,给水度μ=0.03。
图4为解析解、最小二乘配点法和径向基函数配点法三种方法的计算结果。从图中可以看出,随着时间的增加,最小二乘配点法的计算结果比径向基函数配点法更精确,稳定性更好。
图3
图4
上述算例表明:最小二乘配点法处理非承压含水层稳定和非稳定流的水流问题是有效的,既能节省计算量又能保证精度,其近似解和解析解相比的误差比径向基函数配点法更小。此方法同样适用于二维和三维流问题,方法是类似的。
[1]张雄,刘岩.无网格方法[M].北京:清华大学出版社.2004.
[2]王连军.工程地下水计算[M].北京:中国水利水电出版社.2004.
[3]张雄,刘小虎,宋康祖,等.Least-square collocation mesh less method,Ⅰnt.J.Num.Meth.Engng,2001,51(9):1089 –1100.
[4]Kansa E J.Multiquadrics-a scattered data approximation scheme with applications to computational fluid dynamics Ⅰ.Ⅰ.Comp.Math.Applic.1990,19(8/9):147 – 161.
[5]Kansa E J.Multiquadrics-a scattered data approximation scheme with applications to computational fluid dynamics Ⅰ.Comp.Math.Applic.1990,19(8/9):127 – 145.
[6]周德亮,李静.最小二乘配点法解一维非稳定流问题[J].地下水:2012,34(2):20-21.