李小超， 張耀明

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

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

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

文獻(xiàn)[11]提出了一種非線性指數(shù)變換法來處理幾乎奇異積分，并將其應(yīng)用在處理三維彈性問題中，取得了理想的結(jié)果.本文將文獻(xiàn)[11]提出的變換法與三維直接變量規(guī)則化邊界積分方程相結(jié)合，來處理三維位勢問題中的幾乎奇異積分.數(shù)值算例表明，本文算法穩(wěn)定、效率高，即使計(jì)算點(diǎn)非常靠近真實(shí)邊界，使用較少的離散單元就可獲得令人滿意的計(jì)算結(jié)果.

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

考慮三維位勢問題，區(qū)域內(nèi)點(diǎn)y的位勢及其導(dǎo)數(shù)邊界積分方程可寫為

(1)

(2)

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

(3)

式中r為源點(diǎn)y到邊界上場點(diǎn)x的距離.

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

(4)

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

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

(5)

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

距離函數(shù)r2可表示為

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

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

(6)

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

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

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)可直接應(yīng)用高斯求積公式計(jì)算.下節(jié)中的數(shù)值算例驗(yàn)證了本文所提算法的有效性.

4 數(shù)值算例

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

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

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

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

表 1近邊界點(diǎn)位勢u的計(jì)算結(jié)果

表 2近邊界點(diǎn)位勢梯度?u/?x1的計(jì)算結(jié)果

5 結(jié)束語

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

(a)位勢

(b)位勢梯度圖2 400個(gè)內(nèi)點(diǎn)處位勢與位勢梯度的相對誤差曲面

[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ì)算的一種變量替換法[J]. 力學(xué)學(xué)報(bào), 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].力學(xué)學(xué)報(bào)，2013, 45(6): 908-918.