三維位勢直接邊界元法中的幾乎奇異積分

李小超， 張耀明

(山東理工大學 理學院， 山東 淄博 255091)

三維直接邊界元法的有效實施要求準確計算各種核積分.當源點屬于積分單元時，這些核函數是奇異的，對此已有許多有效的處理方法[1-3].另一個重要的問題是源點接近但不在積分單元上的核積分的計算.理論上，這樣的積分是非奇異的，但是它們又表現出類似奇異性的特征，通常不能通過標準的高斯公式準確計算，稱之為幾乎奇異積分.不同于常規積分，幾乎奇異積分需要特別的處理方法.本文著眼于三維邊界元分析中高階幾何單元上的幾乎奇異積分的數值解法研究.

邊界元分析中，一般地認為，準確計算幾乎奇異積分的重要性僅次于奇異積分，近年來引起了許多學者的關注，發展了許多數值處理技術和方法，如剛體位移法[3]、區間分割法[4]、解析與半解析法[5]及各種非線性變換法[6-10]等，文獻[6]對此進行了較全面的綜述.在這些方法中，變換法具有方法簡單、效率高、適用范圍廣等優點，是目前較盛行的方法，包括：三次多項式變換[7]、Sinh變換[8]及指數變換[9-11]等.

文獻[11]提出了一種非線性指數變換法來處理幾乎奇異積分，并將其應用在處理三維彈性問題中，取得了理想的結果.本文將文獻[11]提出的變換法與三維直接變量規則化邊界積分方程相結合，來處理三維位勢問題中的幾乎奇異積分.數值算例表明，本文算法穩定、效率高，即使計算點非?？拷鎸嵾吔?，使用較少的離散單元就可獲得令人滿意的計算結果.

1 二邊界積分方程中的幾乎奇異積分

考慮三維位勢問題，區域內點y的位勢及其導數邊界積分方程可寫為

(1)

(2)

式中q(x)為場點x的法向位勢梯度，u(x)為場點x的位勢.u*(x,y)為三維位勢基本解，q*(x,y)為位勢基本解對邊界外法線的梯度.yk為源點y坐標.離散積分式(1)，(2)所使用的基本核函數有

(3)

式中r為源點y到邊界上場點x的距離.

對于方程(1)與(2)的離散化形式，當計算點y非?？拷e分單元Γe時，計算點y與積分單元Γe之間的距離非常地小，幾乎接近于零，此時方程(1)與(2)中的積分核將產生不同程度的幾乎奇異性，直接應用標準的高斯求積公式不能準確求解.這些幾乎奇異積分可以表示為

(4)

式中r=‖x-y‖2，α>0為常數，f(x,y)是一個規則函數，在不同的積分中可能不同.

2 高階幾何單元的幾乎奇異積分

(5)

這里Nj(ξ1,ξ2)(j=1,…，8)是插值形函數,ξ1,ξ2是無因次坐標

距離函數r2可表示為

r2(ξ1,ξ2)=d2+(ξ1-η1)2g11+

(ξ2-η2)2g22+(ξ1-η1)(ξ2-η2)g12

(6)

3 幾乎奇異積分的變量替換

針對方程(4)中的幾乎奇異積分，考慮如下指數變換

x=d(em1+m2s-1),y=d(en1+n2t-1),

-1≤s≤1,-1≤t≤1

(7)

將式(7)代入式(4)中，可得

(8)

式中

F(s,t)=[1+(em1+m2s-1)2g11(x,y)+

(en1+n2t-1)2·g22(x,y)+

(em1+m2s-1)(en1+n2t-1)g12(x,y)]α

變換后的積分(8)可直接應用高斯求積公式計算.下節中的數值算例驗證了本文所提算法的有效性.

4 數值算例

例1考慮三維球域的熱傳導問題.如圖1所示，該球域半徑為1，表面溫度已知，可以表示為

圖1 單位球域劃分為100個二次曲面單元

該三維球域表面劃分為100個二次曲面單元，邊界量采用二次不連續插值型函數近似.表1、表2分別列出了沿x1方向內點處的位勢u以及位勢梯度?u/?x1的數值結果.可以看出，當計算點距離邊界不是很近時，使用傳統邊界元法(未加變換)或是本文方法均能得到令人滿意的結果；然而，隨著所取得計算點趨于邊界，當計算點到邊界的距離小于0.001時，傳統邊界元法已經失效.另一方面，即使所取的計算點到邊界的距離達到1×10-10，使用本文方法依然可以取得準確的結果.由此表明本文算法在處理幾乎奇異積分是非常有效的.

上，根據坐標θ與φ，均勻選取400個內點進行計算.圖2 (a)、(b)則給出了這400個內點處的相對誤差曲面，位勢u與位勢梯度?u/?x1的平均相對誤差分別為8.7705×10-4、2.2353×10-3.這個數值結果進一步表明本文算法在處理幾乎奇異積分是的準確性與有效性.

表 1近邊界點位勢u的計算結果

表 2近邊界點位勢梯度?u/?x1的計算結果

5 結束語

針對三維位勢直接邊界元法中的幾乎奇異積分，本文利用一種非線性指數變換法，使用高階幾何單元來描述復雜幾何邊界；構造了新的距離函數；利用拓展的變換來消除被積函數的幾乎奇異性.數值算例表明，本文算法穩定、效率高，即使計算點十分靠近邊界，仍可獲得令人滿意的計算結果.

(a)位勢

(b)位勢梯度圖2 400個內點處位勢與位勢梯度的相對誤差曲面

[1] Tanaka M, Sladek V, Sladek J. Regularization techniques applied to boundary element methods[J]. Applied Mechanics Reviews, 1994, 47: 457-499.

[2] Liu Y J. On the simple solution and non-singular nature of the BIE/BEM-a review and some new results[J]. Engng. Anal. Bound. Elem.,2000, 24: 286-292.

[3] Sladek V, Sladek J, Tanaka M. Regularization of hypersingular and nearly singular integrals in the potential theory and elasticity[J]. Int. J. Numer. Meth. Engng, 1993, 36(10): 1 609-1 628.

[4] Gao X W, Yang K, Wang J. An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals[J]. Engng. Anal. Bound. Elem., 2008, 32: 692-696.

[5] Zhang X S, Zhang X X. Exact integrations of two-dimensional high-order discontinuous boundary element of elastostatics problems[J]. Engng. Anal. Bound. Elem., 2003, 28(7): 725-732.

[6] 張耀明, 孫翠蓮,谷巖.邊界積分方程中近奇異積分計算的一種變量替換法[J]. 力學學報, 2008, 40(2): 207-214.

[7] Telles J C F. A self-adaptive coordinate transformation for efficient numerical evaluation of general boundary element integral[J]. Int. J. Numer. Meth. Engng, 1987, 24: 959-973.

[8] Johnston P R, Elliott D. A sinh transformation for evaluating nearly singular boundary element integrals[J]. Int. J. Numer. Meth. Engng.,2005, 62: 564-578.

[9] Zhang Y M, Gu Y, Chen J T.Boundary element analysis of the thermal behaviour in thin-coated cutting tools[J]. Engng. Anal. Bound. Elem.,2010, 34(9): 775-784.

[10] Gao X W, Davies T G. Adaptive algorithm in elasto-plastic boundary element analysis[J]. Journal of the Chinese Institute of Engineers, 2000, 23(3): 349-356.

[11] 張耀明，李小超，Sladek V，等. 三維邊界元法中高階幾何單元上的幾乎奇異積分[J].力學學報，2013, 45(6): 908-918.