謝貴重， 董云橋， 鐘玉東， 周楓林

(1.湖南大學(xué) 機械與運載工程學(xué)院，長沙 410082；2.鄭州輕工業(yè)大學(xué) 機電工程學(xué)院 河南省機械裝備智能制造重點實驗室，鄭州 450002；3.南華大學(xué) 機械工程學(xué)院，衡陽 421001)

1 引 言

在典型的積分類型數(shù)值方法,如邊界元法和擴展有限元法中[1-6]，由于采用奇異基本解或者奇異形函數(shù)，導(dǎo)致被積函數(shù)具有近奇異性。如果不能精確有效地計算這類積分，該類方法的計算精度會受影響甚至會造成錯誤的結(jié)果。

圍繞近奇異積分的精確計算產(chǎn)生了許多方法，如解析和半解析方法[7,8]、自適應(yīng)細分方法[9-11]及非線性變換方法(包括sinh變換和距離變換)[12-17]。解析與半解析方法具有很高精度，然而需要復(fù)雜的數(shù)學(xué)公式推導(dǎo)，限制了該方法的應(yīng)用范圍。自適應(yīng)細分方法是一種穩(wěn)定且通用的方法，但是需要引入額外的積分點，計算量巨大，這種方法在計算超薄型結(jié)構(gòu)時會嚴重降低計算效率。非線性變換也是應(yīng)用較廣的一種方法，在非線性變換中，sinh變換是一種比較有數(shù)學(xué)含義的方法，可以用迭代將積分點集中在近奇異點附近，以達到改善積分精度的目的。

在傳統(tǒng)的sinh變換方法中，首先需要對其進行極坐標變換和坐標變換。文獻[18]指出，在最近點靠近角點或者積分邊界的時候，被積函數(shù)在經(jīng)過極坐標變換后的環(huán)向或經(jīng)過坐標變換后的兩個方向仍然具有很明顯的近奇異性，需要特別關(guān)注。sinh變換結(jié)合環(huán)向變換是一種有效的方法，但是積分單元需要劃分多個三角形，并且在環(huán)向需要2倍的高斯積分點，導(dǎo)致計算效率受到影響。因此，本文采用一種可以將近奇異性在兩個方向進行分離的(α,β)坐標變換法[3]，再根據(jù)復(fù)變函數(shù)極點理論，讓兩個方向的sinh變換分別在極點處實施。最后，給出幾個不同類型單元上近奇異積分的數(shù)值算例，結(jié)果對比表明，雙向sinh變換結(jié)合坐標變換可以明顯提高近奇異積分的計算精度，并且具有較好的穩(wěn)定性，同時讓近奇異積分對最近點的敏感性降低。

2 近奇異積分形式及積分子單元劃分

2.1 近奇異積分形式

在三維邊界元法中，邊界積分形式為

(λ>0)(1)

式中f(F,S)為標量函數(shù)，其數(shù)值積分可以通過標準的高斯積分求得。F為場點，設(shè)其坐標為(x,y,z)。S為源點，設(shè)其坐標為(x0,y0,z0)。λ取值為0.5或1。為了簡單起見，考慮積分區(qū)域為一個落在XY平面上的三角形，設(shè)源點到該三角形的最近點為(x0,y0,0)。將源點到最近點的距離定義為r0。則式(1)可以變形為

(2)

式中由于代入了場點和源點坐標，f(F,S)的形式變換為f1(x-x0,y-y0,z0)。

(3)

式中θ1,θ2和R(θ)可以從文獻[18]找到。

為與文獻[18]的方法進行對比，本文考慮如式(4)類型的積分，

(4)

2.2 與近奇異積分結(jié)合的單元劃分

在近奇異積分單元的劃分中，首先找到源點到積分單元的最近點。如果最近點落在單元內(nèi)部，將四邊形單元劃分為四個積分子單元如圖1(a)所示；如果最近點落在四邊形單元邊上，將四邊形單元劃分為三個積分子單元如圖1(b)所示，故最多分出四個三角形子單元。

3 結(jié)合(α,β)坐標變換的雙向sinh變換法

3.1 (α,β)坐標變換法

從圖2可以看出，一個三角形積分子單元通過(α,β)坐標變換對應(yīng)成了一個標準的四邊形，具體變換如下，

圖2 (α,β)坐標變換系統(tǒng)

(5)

將式(5)代入式(2)得

(6)

式中J=αSΔ為(α,β)坐標變換的雅克比。SΔ為三角形子單元面積的2倍,g(β)=[(x1-x0)+(x2-x1)β]2+[(y1-y0)+(y2-y1)β]2，為了方便，記g(β)=aβ2+bβ+c。式(6)可以進行如下分離，

(7)

3.2 雙向sinh變換及迭代sinh變換法

由式(7)可知，α和β方向的近奇異性都應(yīng)考慮，即需要考慮如下兩種類型的近奇異積分，

(8)

(9)

針對式(8)近奇異積分的處理，可以使用如下sinh變換,

(10)

根據(jù)g(β)的表達式，式(9)可以寫為

(11)

(12)

式中β∈[0,1]?e∈[0,1]

使用此sinh變換的雅可比為

最終得到對式(7)變換之后的積分形式為

(13)

然后對式(13)分別在ξ和e兩個方向采用高斯積分，即可得到最終的數(shù)值積分結(jié)果。

4 數(shù)值算例

將與基于極坐標的sinh變換及其與環(huán)向變換相結(jié)合的結(jié)果進行對比，來驗證本文方法的精確性和有效性，相對誤差定義為

(14)

式中Inum為使用變換的數(shù)值積分方法，Iexa為參考解，可以通過單元細分并使用大量高斯積分點的方法來獲得。接近率定義為

(15)

式中r0為最近距離，Saera為積分單元的面積。

在以下算例中，用符號rsinh表示極坐標下單向sinh變換，rsinh2表示極坐標下單向迭代sinh變換，αsinhβsinh表示(α,β)變換與雙向sinh變換的結(jié)合(即本文方法)，每一塊積分子單元均使用16×16的高斯積分。

4.1 針對線性四邊形單元的近奇異積分

考慮XY空間內(nèi)一個標準的四邊形線性等參元，四個頂點分別為(-1.0,-1.0)，(1.0,-1.0)，(1.0,1.0)和(-1.0,1.0)。最近點在四邊形內(nèi)的參數(shù)坐標分別為(0.9,0.9)，(0.5,0.5)，(0.9,0.0)和(0.98,0.98)，本算例f1=1，計算結(jié)果列入表1和表2。

由表1和表2可知，針對平面單元的近奇異積分，αsinhβsinh 方法要優(yōu)于其他幾種方法，尤其是最近點靠近邊界的情況。

4.2 針對二次四邊形單元的近奇異積分

考慮XYZ空間內(nèi)一個二次的四邊形等參元，其八個頂點分別為(1.0,1.0,0.1)，(-1.0,1.0,0.1),(-1.0,-1.0,0.1)，(1.0,-1.0,0.1)，(0.0,1.0,0.1)，(-1.0,0.0,0.1)，(0.0,-1.0,0.2)，(1.0,0.0,0.2)，(0.0,0.0,0.2)。最近點在四邊形內(nèi)的參數(shù)坐標分別為(0.9,0.9)，(0.5,0.5)，(0.9,0.0)和(0.98,0.98)。本算例f1=1，數(shù)值計算結(jié)果列入表3和表4。

由表3和表4可知，針對二次單元的近奇異積分，相比其他方法，αsinhβsinh方法都可以獲得較好的計算結(jié)果，尤其是最近點靠近邊界時。

表1 λ =0.5，最近點不同時平面單元的計算誤差

表2 λ=1.0，最近點不同時平面單元的計算誤差

表3 λ=0.5，最近點不同時二次單元的計算誤差

表4 λ=1.0，最近點不同時二次單元的計算誤差

5 結(jié) 論

本文提出了一種雙向sinh變換法來計算四邊形單元的近奇異積分。雙向sinh變換方法主要結(jié)合(α,β)坐標變換和sinh變換。首先,利用(α,β)坐標變換來分離α方向與β方向的近奇異性;再基于α方向與β方向的近奇異積分形式，構(gòu)造了sinh變換和迭代sinh變換來消除其近奇異性;最后,通過數(shù)值算例驗證了本文方法的精確性。數(shù)值結(jié)果表明，對于最近點靠近單元邊界時引起的積分子單元頂端大張角和邊長比較大的情況，雙向sinh變換能保持相當高的計算精度和穩(wěn)定性，是一種理想的消除近奇異性的方法。