边界无单元法求解地下水非均质承压稳定流问题

Boundary Element-Free Method to Solve the Problem of Heterogeneous Confined Groundwater Steady Flow

摘要: 将边界无单元法应用于求解地下水非均质承压稳定流问题，给出了用边界无单元法求解地下水非均质承压稳定流问题的方法。将研究区域划分成若干个子区，在每个子区将渗透系数近似表达成常数，从边界积分方程出发，利用改进的移动最小二乘法构造水头近似函数和流量近似函数，在边界节点处生成线性方程组。通过数值算例反映出边界无单元法相较于边界元法具有较高的计算精度。

Abstract: The boundary element-free method is applied to solve the problem of groundwater heterogeneous confined steady flow, and a method for solving the problem of groundwater heterogeneous confined steady flow is given. The study area is divided into several sub-regions, and the permeability coefficient is approximately expressed as the component constant in each sub-region. Starting from the boundary integral equation, the improved moving least squares method is used to construct the water head approximation function and the flow approximation function, and the linear equations are generated at the boundary nodes. The numerical examples show that the boundary element-free method has higher calculation accuracy than the boundary element method.

1. 引言

对于地下水流的研究过程中，为了拟态地下水的流动，早期使用了水均衡法(Water balance method)或水文地质比拟法(Hydrogeological analogy method)。在20世纪初，E.梅勒第一次用解析法(Analytic method)理论原理描述了地下泉水的流动。20世纪50年代，Γ.H.卡门斯基综合了进水量、水流下渗等环境因素，运用有限差分法(Finite difference method)原理描述了地下水动态变化规律。至20世纪80年代，随机模拟法(Stochastic simulation method)作为地下水数值模拟的常用方法在国外逐渐流行起来。

我国有关地下水数值模拟研究始于20世纪70年代中期。1975年，北京大学数学系肖树铁教授等人利用有限元法(Finite element method)原理模拟了河北王窑铁矿的排水问题；同年，山东大学数学系孙讷正教授等人利用有限差分法对山东莱芜铁矿的排水问题进行了数值模拟。肖树铁，孙讷正等人开创了我国地下水数值模拟研究与应用的先河，推动了我国地下水数值模拟相关工作的发展与进步，我国地下水数值模拟的研究自此得到了迅速发展。2003年，程玉民等人通过将改进的移动最小二乘法(Improved moving least square method，简称IMLS法)与边界积分方程法(Boundary integral equation method，简称BIE法)相结合，提出了边界无单元法(Boundary element-free method，简称BEFM) [1] [2] [3] [4] [5]。边界无单元法既是众多无单元方法中的一种，也是一种无网格数值模拟方法。边界无单元法相较于早期Mukherjee等人提出的边界点法，它无需在节点的影响域内布置Gauss点，使求解过程简洁化；与边界单元法相比，边界无单元法省去了离散区域边界这一过程，避免了因在边界上划分单元而带来的前期巨大的工作量，简化了问题的模拟计算过程，节约了问题解决成本，提高了问题求解精度。在边界无单元法中，代入节点的变量值是该点的真实值。边界无单元法具有直观、公式简单、计算结果精度较高、便于引入边界条件等优点，是完全的无网格方法。

本文将边界无单元法应用于求解地下水非均质承压稳定流问题，利用改进的移动最小二乘法构造近似形函数，在边界节点处生成线性方程组。通过数值算例反映出边界无单元法相较于边界元法具有较高的计算精度。

2. 非均质承压稳定流问题的边界积分方程

考虑如下非均质承压稳定流问题：

$\frac{\partial }{\partial x}\left(A\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(B\frac{\partial u}{\partial y}\right)=F,\left(x,y\right)\in \Omega$ (1.1)

边界条件为

$u\left(x,y\right)=\bar{u}\left(x,y\right),\left(x,y\right)\in {\Gamma }_u$ (1.2)

$q\left(x,y\right)=A\frac{\partial u}{\partial x}{n}_x+B\frac{\partial u}{\partial y}{n}_y=\bar{q}\left(x,y\right),\left(x,y\right)\in {\Gamma }_q$ (1.3)

其中： ${\Gamma }_u$ 和 ${\Gamma }_q$ 分别表示本质边界和自然边界，总边界 $\Gamma ={\Gamma }_u+{\Gamma }_q$， $\Omega$ 表示整个研究区域， $F=F\left(x,y\right)$ 表示渗流场 $\Omega$ 内的源，x和y方向的导水系数分别是 $A=A\left(x,y\right)$， $B=B\left(x,y\right)$， $\bar{u}\left(x,y\right)$ 为 ${\Gamma }_u$ 上的水头函数， $\bar{q}\left(x\right)$ ${\Gamma }_q$ 上的流量函数， $\left(n_x,n_y\right)$ 分别表示边界的外法向方向的余弦。

设 $u^{\text{*}}$ 是权函数，对方程(1.1)用加权残量法可得

${\displaystyle \underset{\Omega }{\int }\left[\frac{\partial }{\partial x}\left(A\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(B\frac{\partial u}{\partial y}\right)\right]u^{*}\text{d}\Omega }-{\displaystyle \underset{\Omega }{\int }F{u}^{*}\text{d}\Omega }=0$ (1.4)

对(1.4)式第一项应用格林第二公式可得

${\displaystyle \underset{\Omega }{\int }F{u}^{*}\text{d}\Omega }={\displaystyle \underset{\Omega }{\int }\left(A\frac{{\partial }^2{u}^{*}}{\partial x^2}+B\frac{{\partial }^2{u}^{*}}{\partial y^2}\right)u\text{d}\Omega }+{\displaystyle \underset{\Gamma }{\int }{u}^{*}q\text{d}\Gamma }-{\displaystyle \underset{\Gamma }{\int }{q}^{*}u\text{d}\Gamma }+{\displaystyle \underset{\Omega }{\int }\left(\frac{\partial A}{\partial x}\frac{\partial u^{*}}{\partial x}+\frac{\partial B}{\partial y}\frac{\partial u^{*}}{\partial y}\right)u\text{d}\Omega }$ (1.5)

其中 $u^{*}$ 为方程(1.1)的基本解，即

$\frac{\partial }{\partial x}\left(A\frac{\partial u^{*}}{\partial x}\right)+\frac{\partial }{\partial y}\left(B\frac{\partial u^{*}}{\partial y}\right)+\delta \left(x,y;x_i,y_i\right)=0$ (1.6)

其中 $\delta \left(x,y;x_i,y_i\right)$ 为狄利克雷函数， $\left(x,y\right)$ 和 $\left(x_i,y_i\right)$ 分别表示区域 $\Omega$ 中的动点和定点。

由狄利克雷函数的选择性，式(1.6)可化为

${\displaystyle \underset{\Omega }{\int }F{u}^{*}\text{d}\Omega }=-u_i+{\displaystyle \underset{\Gamma }{\int }{u}^{*}q\text{d}\Gamma }-{\displaystyle \underset{\Gamma }{\int }{q}^{*}u\text{d}\Gamma }+{\displaystyle \underset{\Omega }{\int }\left(\frac{\partial A}{\partial x}\frac{\partial u^{\text{*}}}{\partial x}+\frac{\partial B}{\partial y}\frac{\partial u^{\text{*}}}{\partial y}\right)u\text{d}\Omega }$ (1.7)

其中：令 $x_i=\left(x_i,y_i\right)$，则 $u_i=u(\; x\; i\; )$

对方程(1.5)进行分部积分，可得边界积分方程为

$C_i{u}_i={\displaystyle \underset{\Gamma }{\int }{u}^{*}q\text{d}\Gamma }-{\displaystyle \underset{\Gamma }{\int }u{q}^{*}\text{d}\Gamma }-{\displaystyle \underset{\Omega }{\int }{u}^{*}F\text{d}\Omega }$ (1.8)

其中

$C_i=\{\begin{array}{l}\frac{\theta }{2\pi }当x_i是边界角点时\\ \frac{1}{2}当x_i是光滑边界点时\\ 1当x_i是区域内点时\end{array}$ (1.9)

$u^{*}$ 为(1.1)的基本解， $q^{*}=A\frac{\partial u^{*}}{\partial x}{n}_x+B\frac{\partial u^{*}}{\partial y}{n}_y$

当渗透系数A和B为常量时，由Laplace方程的基本解，可以求出方程(1.1)的基本解 [6] [7]

$\{\begin{array}{l}{u}^{*}=\frac{1}{2\pi \sqrt{AB}}\ln \frac{1}{r_0}\\ r_0=\sqrt{\frac{{\left(x-x_i\right)}^2}{A}+\frac{{\left(y-y_i\right)}^2}{B}}\end{array}$ (1.10)

对于非均质承压稳定流问题，可以将整个区域划分为若干个子区，对于每一个子区，将渗透系数近似表达成分片常数，比如用其形心处的渗透系数代表整个子区的渗透系数。这样，每个子区的基本解由(1.10)式确定。

3. 非均质承压稳定流问题的BEFM法

将研究区域 $\Omega$ 划分为M个子区，在整个研究区域的边界上布置N个节点 $x_i$，通过改进的移动最小二乘法构建近似水头函数和近似流量函数 [8] [9]

$u\left(x\right)\approx {\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_j\left(x\right)u_j}$ (2.1)

$q\left(x\right)\approx {\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_j\left(x\right)q_j}$ (2.2)

其中： $x=\left(x,y\right)$， ${\Phi }_j$ 为改进的移动最小二乘法中的形函数， $u_j$， $q_j$ 分别表示边界节点 $x_i$ 处的水头函数近似值和流量函数近似值。把(2.1)和(2.2)式代入边界积分方程 中确定在每个边界节点处成立

$C_i{u}_i={\displaystyle \underset{\Gamma }{\int }{u}^{*}\left(x_i,x\right){\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_j\left(x\right)q_j\text{d}\Gamma }}-{\displaystyle \underset{\Gamma }{\int }{q}^{*}\left(x_i,x\right){\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_j\left(x\right)u_j\text{d}\Gamma }}-{\displaystyle \underset{\Omega }{\int }{u}^{*}\left(x_i,x\right)F\text{d}\Omega }$ (2.3)

则(2.3)式可写成

$C_i{u}_i+{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left({\displaystyle \underset{\Gamma }{\int }{q}^{*}\left(x_i,x\right){\Phi }_j\left(x\right)\text{d}\Gamma }\right)u_j}={\displaystyle \underset{j=1}{\overset{N}{\sum }}\left({\displaystyle \underset{\Gamma }{\int }{u}^{*}\left(x_i,x\right){\Phi }_j\left(x\right)\text{d}\Gamma }\right)q_j}-{\displaystyle \underset{\Omega }{\int }{u}^{*}\left(x_i,x\right)F\text{d}\Omega }$ (2.4)

根据(2.4)式，将边界上的所有节点 $x_i$ 依次作为定点，此时把所有边界节点得到的方程联立，可列出代数方程组

$\left[C_{\Gamma }\right]\left[U_{\Gamma }\right]+\left[\tilde{H}\right]\left[U_{\Gamma }\right]=\left[G\right]\left[Q_{\Gamma }\right]+F$ (2.5)

其中

$\begin{array}{l}{U}_{\Gamma }={\left(u_1,u_2,\cdots ,u_N\right)}^{\text{T}},Q_{\Gamma }={\left(q_1,q_2,\cdots ,q_N\right)}^{\text{T}},\\ {\tilde{H}}_{ij}={\displaystyle \underset{{\Gamma }_j}{\int }{q}^{*}\left(x_i,x\right)}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_J^{*}\left(x\right)\text{d}\Gamma },G_{ij}={\displaystyle \underset{{\Gamma }_j}{\int }{u}^{*}\left(x_i,x\right){\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_J^{*}\left(x\right)\text{d}\Gamma }},\\ F_i=-{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\displaystyle \underset{{\Omega }_j}{\int }{u}^{*}\left(x_i,x\right)F\text{d}\Omega }},C_{\Gamma }=\left(0,\cdots ,0,C_i,0,\cdots ,0\right)\end{array}$ (2.6)

令

$\left[H\right]+\left[C_{\Gamma }\right]+\left[\tilde{H}\right]$ (2.7)

则(2.7)式变形为

$\left[H\right]\left[U_{\Gamma }\right]=\left[G\right]\left[Q_{\Gamma }\right]+F$ (2.8)

其中：H和G分别是边界点函数值 $U_{\Gamma }$ 和边界点法向导数值 $Q_{\Gamma }$ 的系数矩阵，且H和G都是N维矩阵。 $U_{\Gamma }$ 、 $Q_{\Gamma }$ 及F是N维列向量。

根据边界条件可知，在本质边界上的边界点，u已知而q未知，在自然边界上的边界点，u未知而q已知，所以，在边界上节点 $x_i$ 处有N个未知量对应N个线性方程。因而(2.8)式是含有N个未知量的线性方程组，可矩阵化表示为

$\left[H_1{H}_2\right]\left(\begin{array}{c}{U}_1\\ U_2\end{array}\right)=\left[G_1{G}_2\right]\left(\begin{array}{c}{Q}_1\\ Q_2\end{array}\right)+F$ (2.9)

将(2.9)式转化为

$\left[-G_1{H}_2\right]\left(\begin{array}{c}{Q}_1\\ U_2\end{array}\right)=\left[-H_1{G}_2\right]\left(\begin{array}{c}{U}_1\\ Q_2\end{array}\right)+F$ (2.10)

其中： $Q_1$ 和 $U_2$ 为待解量

在(2.10)式中，用N维的可逆矩阵A表示等式左端的系数矩阵，N维的列向量X表示等式左端的未量，等式右端用N维列向量B表示，可得到线性方程组为

$AX=B$ (2.11)

这样，通过上面的线性方程组(2.11)就能够解出边界节点对应的未知的函数值和未知法向导数值。

当边界节点解的近似值得出后，区域内任意一点 $x_i$ 处解的近似值由(1.4)得到

$u_i={\displaystyle \underset{\Gamma }{\int }{u}^{*}\left(x_i,x\right){\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_j\left(x\right)q_j\text{d}\Gamma }}-{\displaystyle \underset{\Gamma }{\int }{q}^{*}\left(x_i,x\right){\displaystyle \underset{j=1}{\overset{N}{\sum }}{\Phi }_j\left(x\right)u_j\text{d}\Gamma }}-{\displaystyle \underset{\Omega }{\int }{u}^{*}\left(x_i,x\right)F\text{d}\Omega }$ (2.12)

4. 数值算例

考虑一个满足方程 $\frac{\partial }{\partial x}\left(A\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(B\frac{\partial u}{\partial y}\right)=F,\left(x,y\right)\in \Omega$ 的研究区域，这里取 $\Omega$ 为： $0.6<x,y<1.4$，

$M=16$， $N=16$，给出三例：

(1) $A=x$， $B=y$， $F=2$ ；

(2) $A=x+y$， $B=x+y$， $F=2\cos \left(x+y\right)-2\left(x+y\right)\sin \left(x+y\right)$ ；

(3) $A=x$， $B=y$， $F=\cos x-x\sin x$ ；

根据表1~3的计算结果表明该方法可靠，且相较于边界元法具有较高的计算精度。

Table 1. The u value at the centroid of the numerical example (1) is listed and compared with the exact solution and the BEM solution

表1. 列出数值算例(1)子区内形心处u值并与准确解和边界元法解作比较

Table 2. The u value at the centroid of the numerical example (2) is listed and compared with the exact solution and the BEM solution

表2. 列出数值算例(2)子区内形心处u值并与准确解和边界元法解作比较

Table 3. The u value at the centroid of the numerical example (3) is listed and compared with the exact solution and the BEM solution

表3. 列出数值算例(3)子区内形心处u值并与准确解和边界元法解作比较

5. 结语

本文提出了用边界无单元法求解非均质承压稳定问题的方法，通过数值算例表明该方法计算结果十分理想。将边界无单元法与边界元法进行了比较，精度较边界元法有所提高。该方法还具有节点布置灵活，完全不需要网格等优点，但由于边界无单元法需要进行大量的积分运算，这可能会使计算结果产生误差。如何改进该方法，将有待进一步研究。