基于坐標變換精確計算近奇異積分的雙向sinh變換法

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

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

1 引 言

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

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

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

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

2.1 近奇異積分形式

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

(λ>0)(1)

式中f(F,S)為標量函數，其數值積分可以通過標準的高斯積分求得。F為場點，設其坐標為(x,y,z)。S為源點，設其坐標為(x0,y0,z0)。λ取值為0.5或1。為了簡單起見，考慮積分區域為一個落在XY平面上的三角形，設源點到該三角形的最近點為(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 與近奇異積分結合的單元劃分

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

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

3.1 (α,β)坐標變換法

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

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

(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)可知，α和β方向的近奇異性都應考慮，即需要考慮如下兩種類型的近奇異積分，

(8)

(9)

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

(10)

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

(11)

(12)

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

使用此sinh變換的雅可比為

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

(13)

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

4 數值算例

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

(14)

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

(15)

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

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

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

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

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

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

考慮XYZ空間內一個二次的四邊形等參元，其八個頂點分別為(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)。最近點在四邊形內的參數坐標分別為(0.9,0.9)，(0.5,0.5)，(0.9,0.0)和(0.98,0.98)。本算例f1=1，數值計算結果列入表3和表4。

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

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

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

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

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

5 結 論

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