在边界单元法中,由于积分核的奇异性,采用通常高斯求积公式是无法精确计算靠近边界的场量,这一现象即所谓边界层效应.避免边界层效应的间接方法有刚体位移法[1]、补域法[2]和插值法[3]等,直接方法有子单元法[1]等.刚体位移法对于裂纹问题并不适用[4],也不适于场变量的导函数的计算.补域法需要延拓所关心区域附近的边界以便形成补域和虚拟的边界,而插值法则需选用适当的插值函数并计算边界点与内点的场变量.子单元法的不足之处在于当计算点离边界的距离较小时,所需剖分的子单元数目较大,采用低次单元时的效果并不理想.
本文以二维位势问题为例,给出了解决边界单元法中的边界层效应并精确计算边界近旁场量的一般方法,这种方法与边界的离散方式无关.
1 方法介绍
1.1 位势问题的积分方程
对于以Γ为边界的区域Ω,描述位势问题的场量及其导函数的积分方程如下:
(1)
(q∈Ω)(2)
式中 u位势
n边界外向法线
p,q源点与场点
u二维位势问题的基本解
(3)
式中,r表示p,q两点间的距离.当p点距离边界Γ较近时,由于积分核表现出奇异性,位势u及其导函数u,k是不能直接用式(1)和式(2)通过高斯求积公式准确计算的,此即所谓边界层效应.其中式(1)右边第一个积分核具有对数奇异性,式(1)右边第二个积分核与式(2)右边第一个积分核具有Cauchy奇异性,式(2)右边第二个积分核具有Hardama奇异性.
1.2 距离函数的定义
准确计算奇异积分的关键是在单元的自然坐标系中把积分核的奇异部分分离出来[4].设边界Γ已被剖分成若干个边界单元,Φ(ξ)为单元坐标变换的某一形函数,G(ξ)为变换Jacobi,ξ为局部坐标,对于靠近p点的单元积分的一般表达式为
(4)
式中,KS(p,q)表示某一积分核,下标S代表积分核的种类,其中S=L代表对数核,S=C代表Cauchy核,S=H代表Hardama核.首先定义r0为p点到单元的最小距离,ξ0为单元上对应于最小距离的点的局部坐标.有两种情况需要区分,如图1所示,对于单元AB,最小距离为pq0,ξ0对应于q0点.而对于单元BC,最小距离为pB,这时ξ0对应于单元的端点B.
图1 p点近旁的边界单元
对于靠近p点单元的距离函数可以定义为
(5)
1.3 奇异部分的分离
借助于距离函数,便可以分离出积分核的奇异部分:
(6)
式中,FS[g(ξ)]分别定义为
(7)
(8)
(9)
以上3式分别对应于具有对数奇异、Cauchy奇异和Hardama奇异性积分核的奇异部分.式(6)中右边的第一个积分是正常积分,可用高斯求积公式准确计算.右边的第二个积分可积,
即
(10)
(11)
(12)
以上式中:
G0=G(ξ0) (13)
这样,根据式(6)、式(10~12),采用通常的高斯求积公式就可以准确计算靠近p点单元的积分值,有效地解决边界单元法中存在的边界层效应问题.
2 算例与讨论
采用三次Hermite样条插值函数,即三次单元[5]对椭圆区域上的位势问题进行了计算,椭圆长半轴a=2,短半轴b=1.椭圆的参数方程为
x1=asinθ (14)
x2=bcosθ (15)
在椭圆上划分了16个边界单元,每个象限4个边界单元.采用的位势定义为
u=ln(R) (16)
其中:
(17)
在计算中取x01=-1,x02=3.对普通单元采用4点高斯积分,对p点近旁单元采用8点高斯积分.计算结果列于表1~4,表中,p(xp1,xp2)表示p点的位置坐标,r0表示p点到边界的最小距离,r0/l为相对距离,其中l为单元长度.计算实践表明,在计算位势的导函数值时,由于Hardama积分核的超强奇异性(式(2)右边第二个积分),对于如图1所示的单元AB的情况,仍然需要但仅仅需要在q0点剖分一次子单元,才能保证必要的计算精度.
表1 计算值与理论值的比较(p点位于第一象限)