科学与工程领域中许多问题,例如,稳定热传导、弹性杆件的扭转、稳定渗流、动水压力、薄膜平衡、Helmhotz方程、电磁场、非均质材料及非线性问题等^[1-3],的数值分析可直接地或者间接地归结为Laplace或Poisson方程的边值问题的求解. 边界元法是求解此类问题的强有力的数值工具.

物理量梯度边界积分方程由物理量边界积分方程关于空间变量取导数获得^[4],它对于某些问题,如裂纹问题、波散射问题、薄板弯曲问题及具有退化边界的狭窄和薄域问题等^[5-7],是非常有用的,因为此时仅由基本物理量的边界积分方程无法准确地表示原边值问题的解,即与原边值问题不等价.为了避免出现这种情况,通常将基本物理量的梯度边界积分方程与物理量边界积分方程组合,即对偶边界积分方程,来表示原边值问题的解.位势问题的位势梯度边界积分方程已收到了许多研究^[8-12].文献^[8]建立了二维位势问题的直接变量规则化梯度边界积分方程,文献^[9]建立了三维位势问题的直接变量规则化梯度边界积分方程.两个梯度方程中,超奇异积分的规则化形式是相似的,即

$ \int_\varGamma {\left( {u({\pmb x}) - u({\pmb y}) - u_{,k} ({\pmb y})(x_k - y_k } \right)O \Big (\dfrac{1}{r^2} \Big ) d\varGamma _{\pmb x} } $

需指出的是,$u_{,k} ({\pmb y})$既不是已知量,也非未知量,数值实施时需将$u_{,k} ({\pmb y})$沿曲面(曲线)的两个(一个)切线方向及法线方向进行分解,切线方向的量采用插值逼近,法线方向的量作为未知量,因此这个积分在高价曲面单元下的数值实施是非常复杂和困难的.文献^[10]研究了二维位势问题的位势梯度边界积分方程.引入一系列变换将梯度边界积分方程转换为具有强奇异积分的自然边界积分方程,然后采用加减技术规则化强奇异积分.方法有效地消去了奇异性,取得了很好的数值结果,但很难推广到三维位势问题.文献^[11,12]建立了三维位势问题的间接变量规则化梯度边界积分方程,给出了高阶单元下的数值实施方案,取得了准确的数值结果.方程中规则化积分形式是

$ \int_\varGamma \Big [ u({\pmb x}) - u({\pmb y}) \Big ] O \Big (\dfrac{1}{r } \Big ) d\varGamma _{\pmb x} $

它比直接法中的形式简单得多,数值实施也容易得多.但是,其数值实施也需要许多技术性措施^[11-14].不同于前述的规则化边界积分方程法,文献^[15,16,17]和文献^{[18,19,20,21]}分别提出了直接计算梯度方程中的奇异积分的局部规则化方法.它们的优点是算法具有一般性,可计算任何阶的奇异积分. 缺点是两种算法的理论推导是繁复的,计算量也相当大,因为它们需将积分核中的每个依赖于单元参数的函数表示成在内蕴坐标 系统下的局部距离$\rho $的幂级数.此外,还有一些别的直接计算方法^[22-26],都属于局部法,涉及与局部单元参数有关的操作.

本文提出了求解三维位势问题的位势梯度边界积分方程的新算法.该方法通过构造辅助边值问题来计算梯度边界积分方程中的系统矩阵,算法实施中不需要建立规则化边界积分方程也不需要直接计算强奇异积分.因此,方法具有数学理论简单、计算效率高、结果精度高等优点. 需要强调的是,本文辅助边值问题法为梯度边界积分方程中的强奇异积分计算提供了一个新的思路和途径.另外,本文方法可以拓展到其他问题,如弹性力学问题、Stokes及Helmholtz问题等.

本文设$\varOmega$是示$R^3$中的一个有界区域,$\varOmega^c$是它的补域,$\varGamma$是它 们的边界.${\pmb n}=(n_1,n_2)$是区域$\varOmega$的边界$\varGamma$在$x$点处的单位外 法向量.三维位势问题的基本解是

区域$\hat\varOmega \subset R^3(\hat\varOmega =\varOmega$ 或$\varOmega^c$)上的三维位势问题的控制微分方程(不含源项)是

具有边界条件

这里$\varGamma =\varGamma_1 \cup \varGamma_2$且$\varGamma_1 \cap \varGamma_2=\phi$; $\bar u({\pmb x})$和$\bar q({\pmb x})$是边界$\varGamma$上已知的函数.

引理^[27-30] 设$\psi({\pmb x})\in C^{0,\alpha} (\varGamma)$和$\hat{\pmb x}$是曲面$\varGamma$上的任一光滑点,当${\pmb y} \to \hat{\pmb x} \in \varGamma$时,令$h=|{\pmb y}-\hat{\pmb x}|$,$d=\mathop{\rm inf}\limits_{x \in \varGamma} | {\pmb y}-{\pmb x}|$和假设$\dfrac hd\leq K_1$ ($K_1$是一个常数),那么有

区域$\hat{\varOmega} \subset R^3$内的位势函数可表示为

这里$\hat{\varOmega}=\varOmega$或$\varOmega^c$. 令${\pmb y} \to \varGamma $,可获得如下边界积分方程

现在从式(6)出发,导出位势梯度边界积分方程. 对任一个给定点$\hat {\pmb x} \subset \varGamma$,由式(6)可得

让${\pmb y} \to \hat{\pmb x}$,利用定理1,可得如下位势梯度边界积分方程(将$\hat{\pmb x}$用${\pmb y}$代替)

这里$\begin{cases} S=\left\{ 1 , \ \hat{\varOmega}={\varOmega} \\ -1 , \hat{\varOmega}={\varOmega}^c\!\!\right.\end{cases}$,C.P.V.表示柯西主值积分.

边界元法的实施包括边界几何的描述和边界函数的插值逼近. 为了一般性,本文采用八节点四边形二次单元来描述边界几何,八节点二次不连续插值函数来逼近单元上的边界函数.

八节点四边形二次单元的形函数$N^k(\xi_1,\xi_2)$ $(k=1,2,\cdots,8)$是

这里$\xi_1,\xi_2$是无因次坐标,且$-1\leq \xi_1$, $\xi_2 \leq 1$. 因此单元上的点${\pmb x}$可以近似地表示为

这里${\pmb x}^k$是第$k$个插值节点的坐标.

单元上的边界量由八节点二次不连续插值函数来逼近,即

这里$\phi^k$是边界函数在第$k$个节点处的函数值,$\alpha \in (0,1)$ 的常数,本文取参数$\alpha=2/3$.

许多情况下,边界几何具有参数表示. 此时采用精确单元计算,可减少误差. 精确单元的概念最早是本文作者2004年提出的^[27],后来许多学者跟随了这个方法,甚至赋予一个新的名字,但本质是一样的.

在三维空间中,多数情况下,固体模型的表面可以表示成参数形式

这里$a, b, c, d$都是有限数. 当区间$[a,b]$ 和$[c,d]$ 分别被离散成$N_1$和$M_1$个子区间后,相应的曲面被离散成$ N\times M$个所谓的单元. 由于这样的分割是在参数空间内,对应的几何点仍然在原来的曲面上,因此称为精确单元. 在每个精确单元上,参数$\theta,\varphi$可表示为

这里$\theta_i,\varphi_i \ (i=1,2)$ 是精确单元的端点坐标.

将边界$\varGamma$离散成$N_{\rm e}$单元,因而有$8\times N_{\rm e}$个边界节点：${\pmb x}^i$,$i=1,2,\cdots,8\times N_{\rm e} $. 方程(7)和(9)中的${\pmb y}$取为任一场点${\pmb x}_i \in \varGamma_{\alpha}$(属于第$\alpha$个单元),那么方程(7)和(9)的离散化形式是

这里

$ G_{jl}=\int_{\varGamma_j} N^l ({\pmb x}) u^* ({\pmb x},{\pmb x}^i) d \varGamma = \\ \qquad \int^1_{-1} \int^1_{-1} N^l ({\pmb x}(\xi_1,\xi_2)) u^*({\pmb x}(\xi_1,\xi_2),{\pmb x}^i) J(\xi_1,\xi_2) d \xi_1 d \xi_2\\ Q_{jl}=\int_{\varGamma_j} N^l ({\pmb x}) \dfrac{\partial u^*({\pmb x},{\pmb x}^i)}{\partial y_k} d \varGamma= \\ \qquad \int^1_{-1} \int^1_{-1} N^l ({\pmb x}(\xi_1,\xi_2)) \dfrac{\partial u^*({\pmb x}(\xi_1,\xi_2),x^i)}{\partial y_k} J(\xi_1,\xi_2) d \xi_1 d \xi_2 $

其中, $ J(\xi_1,\xi_2)=|{\pmb P}_1\times {\pmb P}_2 |$,且

$ {\pmb P}_1=\Big ( \dfrac{\partial x_1}{\partial \xi_1}, \dfrac{\partial x_2}{\partial \xi_1}, \dfrac{\partial x_3}{\partial \xi_1} \Big ) , \ \ {\pmb P}_2=\Big ( \dfrac{\partial x_1}{\partial \xi_2}, \dfrac{\partial x_2}{\partial \xi_2}, \dfrac{\partial x_3}{\partial \xi_2} \Big ) $

数值实施中,确定式(16)中的$A^i$是一个困难的问题.$A^i$对应离散系统矩阵中的对角元,并且对角占优,因此$A^i$的准确与否对解系统的性能及解的精度影响特别大.另一方面,估计$A^i$需要计算一个强奇异积分,其强奇异核函数是基本解关于坐标变量的偏导数,一般不是边界法向导数,因此它的计算是相当困难的^[11-12].采用规则化边界积分方程可避免直接计算此类积分^[11-14],但规则方程的形式较复杂,程序设计难度较大;直接计算此类积分需要繁复的理论推导和计算工作量^[15-21]. 本文提出计算式(15),式(16)对应的系统矩阵$[{ G}_{jl}]$和$[{ Q}_{jl}]$的辅助边值问题法.构造辅助超定边值问题,求解此边值问题可求得$[{ G}_{jl}]$和$[{ Q}_{jl}]$. 下面给出具体实施过程:

(1) 当$\hat \varOmega=\varOmega$时,构造如下超定辅助边值问题

这里$\bar{u}$和$\bar{q}$由下面的函数确定

(2) 当$\hat \varOmega=\varOmega^c$时,构造如下超定辅助边值问题

这里, $\bar{u}$和$\bar{q}$由下面的函数确定

这里$(a,b,c)$ 是$\varOmega^c$外的任一点.

现在用式(15)和式(16)求解边值问题(17)或(19),就可获得$[{ G}_{jl}]$和$[{Q}_{jl}]$.值得注意的是,在求出这两个矩阵后,使用式(15)和式(16)求解有限域$\varOmega$或者无限域$\varOmega^c$上的任何边值问题,不再需要计算$[{ G}_{jl}]$和$[{ Q}_{jl}]$. 由此可看出,辅助边值问题法的效率并不差.

考虑3个数值算例,来验证方法的有效性. 数值实验的重点在于考察辅助边值问题法计算边界通量${\partial u} /{\partial x_1 }$的能力及准确性. 为了估计数值误差,采用如下$L_2 $范数

其中, $M$表示计算点数,$u_{\rm exact}^k $和$u_{\rm num}^k $分别是第$k$个计算点处的精确解和数值解.$RE_u $,$RE_q $及$RE_{q_1 } $分别表示边界位势$u$、法向梯度$q = \partial u / \partial {\pmb n}$及热流通量$q_1 =\partial u / \partial x_1 $的平均相对误差.

算例 作为第一个算例,考虑立方体区域$\varOmega=\{ (x_1,x_2,x_3) \in R^3: 0 \leq x_1,x_2,x_3 \leq 1 \}$上的热传导问题,如图1所示. 边界条件如下

$ \left\{\begin{array} u\left( {x_1 = 0} \right) = {x_2^2 } / 2 - x_3^2 - 5x_2 - 6x_3 \\ u\left( {x_1 = 1} \right) = {x_2^2 } / 2 - x_3^2 - 5x_2 - 6x_3 - 7 / 2 \\ u\left( {x_2 = 0} \right) = {x_1^2 } / 2 - x_3^2 - 4x_1 - 6x_3 \\ u\left( {x_2 = 1} \right) = {x_1^2 }/ 2 - x_3^2 - 4x_1 - 6x_3 - 9 / 2 \\ q\left( {x_3 = 0} \right) = 6,\quad q\left( {x_3 = 1} \right) = - 8 \end{array}\right. $

这是文献^[11]采用的算例,文中没有给出边界量的计算结果,因而这里不便于比较. 立方体的边界划分为54个四边形线性单元,每一正方形表面划分9个单 元,如图1所示. 图2给出了$x_3 =1$面上的直线$x_1 = {11}/{18}$上的9个边界节点处的 温度$u$和热通量$q_1 = {\partial u}/{\partial x_1 }$的数值计算结果,可看出,数值解和精确解 吻合的很好.此外,在立方体内部选取400个均匀分布在区域 $\left\{ {0.15 \leqslant x_1 ,x_2 \leqslant 0.85,x_3 = 0.5}\right\} $上的计算点,图3(a)和图3(b)分别展示了400个计算点上的温度$u$和通 量$q_1 = {\partial u} / {\partial x_1 }$数值解的相对误差曲面.

算例2 为了与文献^[11]提出的规则化边界元法进行比较,本算例采用文献^[11]中的算例. 考虑单位球上的热传导,如图4所示. 问题具有Dirichlet边界条件：$\left. {u(r,\varphi ,\theta )} \right|_\Gamma = \bar{u}$,其中$\bar {u}$由下列问题的解析解给出

$ u = {\sin ^2\varphi \left( {2\cos ^2\theta - 1} \right)} /2 + 6\sin \theta \cos \varphi + 8\cos \varphi $

在单位球内部选择400个计算点,均匀分布在球面(根据$\theta ,\varphi $) $\left\{ {\left( {x_1 ,x_2 ,x_3} \right) \in R^3:x_1^2 + x_2^2 + x_3^2 = 0.5} \right\}$上.图5描述了这400个计算点上的场温度$u$和场梯度${\partial u} / {\partial x_1}$数值解的$L_2$误差随边界单元数增加的变化情况,即收敛曲线.表1列出了本文方法与规则化边界积分方程法^[11]在不同单元数下,求得的边界位势梯度$q = {\partial u} /{\partial x_1 }$和法向梯度$q = {\partial u} / \partial {\pmb n}$数值解的相对误差$RE_{q_1 }$和$RE_q $,从对比中可看出,辅助边值问题法比规则化边界积分 方程的准确性稍好,这主要是因为辅助边值问题法计算系数矩阵的对角元可能更准确.

算例3 在最后一个算例中,考虑一个复杂区域上的热传导问题. 计算区域是如图6所示的一个变形长方体^[31],它是由长方体$ \{ L\times W\times H =$ $5.0\times 1.0\times1.0\}$变形得到,其中$L$,$W$,$H$分别代表长方体的长、宽、高. 变形长方体的上、下面的两条长边可分别表示为$x_3 = 0.1x_2^2 $和$x_3 = 0.1x_2^2 + 1.0$,左、右侧面的宽度和高度不变.在变形长方体边界上施加混合边界条件,其中变形长方体下面通量$\bar {q}$已 知,其余各面温度$\bar{u}$已知$\bar {u} = x_1 x_2 x_3 + 10x_1 + 10x_2 + 10x_3 $, $\bar {q} = (x_2 x_3 + 10)n_1 + (x_1 x_3 +10)n_2 + (x_1 x_2 + 10)n_3 $问题的解析 解为$u(x_1 ,x_2 ,x_3 ) = x_1 x_2 x_3 + 10x_1 + 10x_2 + 10x_3 $.

采用混合单元计算. 变形长方体的边界被划分为48个单元,其中左、右侧面各划分为4个线性单元(精确单元),其余4个面各划分为10个8节点四边形2次单元,如图5所示. 沿着变形长方体的中心线$\left\{ {x_3 = 0.1x_2^2 + 0.5,x_1 = 0.5}\right\}$选取10个计算点,表2列出了10个计算点上的温度$u$数值解以及精确解.此外,图7给出了 变形长方体的5个表面(不包括下表面)上的边界节点处的法向梯度${\partial u} /{\partial {\pmb n}}$和 梯度${\partial u} /{\partial x_1 }$的相对误差随边界单元数增大的变化情况,即收敛曲线.

本文提出求解位势梯度边界积分方程的辅助边值问题法. 构造与边值问题具有相同解域的辅助边值问题,通过求解辅助边值问题,可准确地获得梯度边界积分方程的离散系统矩阵,计算过程不涉及强奇异边界积分的计算. 所求得的系统矩阵可直接用来求解边值问题,不再需要重新计算系统矩阵. 三个数值算例验证了方法的可行性、准确性及收敛性. 此外,从与规则化边界元法的数值结果的对比可看出,本文方法的数值结果稍好一点,说明采用辅助边值问题法计算梯度边界积分方程的系统矩阵更准确,特别是对角元的计算.

本文为梯度边界积分方程的求解提供了新的思路. 辅助边值问题方法具有理论简单,程序设计容易,精度高等优点,而且容易拓展到弹性问题、Stokes问题、Helmholtz问题等.