精確計算彈性薄型結構中奇異積分方法

郭壯志

（1.鄭州輕工業大學機電工程學院,河南 鄭州 450002；2.河南省機械裝備智能制造重點實驗室,河南 鄭州 450002）

在彈性問題的邊界元法求解中，弱奇異積分一直是其研究的一個重要課題[1-4]。尤其對于薄型結構[5-8]，此結構在進行離散時，將會存在形狀較差的單元（比如含大角度或狹長的單元），嚴重影響計算結果的精度。準確快速的計算邊界積分方程中的弱奇異積分是其順利求解的關鍵。

目前有很多處理弱奇異積分的方案，比如積分簡化法、單元細分法和坐標變換法等[9-12]，每一種方法都存在優缺點。但當源點位置比較差時（比如源點靠近單元的邊界），單純的積分變換無法獲得較高的計算精度，需要對單元進行細分。因此，單元細分法加上極坐標變換方法是最廣泛應用的一種方法。目前存在不少的細分方法，比如奇異點（即源點）直接和單元頂點相連、區間分塊法，等等，這些方法在進行細分的時候，會分出一些形狀比較差的積分塊（比如含大角度的三角形積分塊和邊長比較大的積分塊），為奇異積分的處理增加了困難。極坐標變換法是消除弱奇異性的關鍵，它把面積分轉換為環向和徑向的雙重積分，再對徑向作變換來達到消除奇異性的目的。但這種方法在使用時，每次都需要重新計算積分塊在徑向和環向的積分區間，實現起來比較的麻煩。

本文針對彈性薄型問題邊界積分方程中位移基本解的弱奇異性質，對于一個離散單元的數值積分，首先依據源點位置、單元形狀和源點到單元的最近距離，分別開發三角形和四邊形單元的細分技術；然后在該細分技術的基礎上，針對細分所得到的積分塊，構造一種更加簡單的坐標變換法，用來消除被積函數中的奇異性。該方法相比常規的極坐標變換法，不需要再去計算它們的積分區間，實現起來更加簡單有效；最后將該細分技術和坐標變換方法用于彈性薄板結構的邊界元法求解，實現該問題精確有效的求解。

1 薄型問題的奇異邊界積分方程

三維彈性薄型問題的邊界積分方程如式（1）所示，其中P和Q分別為源點和場點，在光滑邊界上cij(P)=0.5，和為位移和面力基本解，其表達式如式（2）所示。

其中，G和v分別是剪切模量和泊松比，r是源點和場點之間的距離，uj(Q)和tj(Q)分別表示邊界上的位移和面力，δij為克羅內克符號，ni和nj為i、j方向的外法向量。若把方程(1)中的邊界離散成N個單元，式（1）可以寫為：

其中，n是每一個單元上的節點數，Nα(Q)是單元上第α個節點的形函數。從式(2)可以看出，對位移基本解的積分具有弱奇異性，而對面力基本解的積分具有強奇異性，本文重點關注對位移基本解的積分。

2 弱奇異邊界積分處理方法

首先取一個離散單元作為研究對象，把對位移基本解的積分簡化成式（4）形式，其中，P和Q分別代表著源點和場點，r是P和Q之間的距離，f(P,Q)為是一個不帶奇異性的光滑函數，Ф(Q)是相應的形函數。

針對此離散單元的積分具有弱奇異性，下面介紹一下對該單元的奇異積分處理方法，以達到精確求解弱奇異積分的目的。先來介紹一種奇異積分單元細分技術。

2.1 求解奇異積分的單元細分技術

由于源點位置的不確定性（可能落在單元的內部、邊和頂點上），所以細分方式也會不同。當源點位置給定時，首先計算源點到單元各邊的最近距離d，然后以0.2 d為長度在源點附近構造一個很小的四邊形，并把源點包括在里面。再分別延長小四邊形的四條邊，分別與積分單元進行相交，得到若干個三角形和四邊形積分塊。最后把包含源點的小四邊形的各個頂點與源點相連得到若干個含奇異點s的三角形積分塊。三角形和四邊形單元的具體細分方案如圖1和圖2所示。其中，0.2 d的選取是為了避免積分塊之間相互干涉，根據數值試驗得到的一個經驗數值。

圖1 不同源點位置三角形單元的細分方法

圖2 不同源點位置四邊形單元的細分方法

除此之外，若積分單元比較狹長，可以先根據該單元的最長(Lmax)和最短邊(Lmin)的比例先對其進行一次細分（如圖3所示，若Lmax/Lmin＜ 2則不分；若2 ＜Lmax/Lmin＜ 3，就把該單元分成兩個積分塊。以此類推），然后按照圖1和2中的細分方式對含源點的積分塊進一步細分。另外，若細分得到細長的積分塊（圖1和2中細長的四邊形積分塊），仍可根據圖3中的方法，對該積分塊進一步細分，以確保細分結果中不出現含較大邊長比的積分塊。

圖3 細長四邊形單元的細分方法

2.2 弱奇異性消除方法

通過上述單元細分可知，不管是三角形還是四邊形單元，最后都會把包含奇異點（源點）的小四邊形分成若干個三角形，而源點剛好是這些三角形的其中一個頂點。不含源點的積分塊可以通過普通的高斯積分（四邊形積分塊）或漢默積分（三角形積分塊）來求解，但由于包含源點的三角形積分塊具有奇異性，單純的漢默積分無法消除被積函數中的奇異性，需要進行特殊處理。

以式(4)作為積分方程，針對細分得到的包含源點的積分塊，本文構造一種新型的坐標變換法來消除其奇異性。引入了(α,β)局部坐標系，相應的坐標變換如圖4所示。

圖4 三角形積分塊的坐標變換

坐標變換的具體過程如下：

將式（5）和（6）代入（7）可得

通過式（8）的坐標變換，把原自然坐標系下的三角形映射為局部坐標系下的一個正四邊形，四邊形的兩個頂點坐標分別為（0，0）和（1，1），即局部坐標α和β的取值范圍為[0,1]。坐標變換的雅可比為αS。

經過變換之后，式（4）在奇異三角形積分塊的積分就可以寫成如下形式：

從式（5）～式（9）可以看出，這種新型的坐標變換得到的雅可比中含有擬零因子α，剛好和積分方程（4）分母中的擬零因子抵消掉，達到消除奇異性的目的。并且這種坐標變換方法直接把α和β的積分區間變換到[0,1]，不需要再去計算每一個積分塊的積分區間，相比傳統的極坐標變換實現起來更加的簡單有效。

3 數值算例

本節給出了幾個算例來驗證本文方法的精確性和有效性。首先驗證本文方法在一個單元上的積分精度（算例考慮了不同源點位置和單元的形狀類型），最后把本文方法用于彈性薄板結構的求解。相對誤差的計算式如式(12)所示。

其中，In為數值計算結果，Ie為相應的參考解。參考解通過使用單元細分和坐標變換法，采用大量高斯積分點得到。為了方便計算，假設在式(4)中，f(P,Q)=1.0。下述第一個例子考慮式（4）形式具有弱奇異性的邊界積分。

3.1 奇異積分算例

本例以四邊形單元為例來驗證本文方法。四邊形單元的四個頂點為(0.0,0.0,0.0)、(1.0,0.0,0.0)、(1.0,1.0,0.0)和(0.0,1.0,0.0)，邊長比為a=1。采用本文方法和傳統方法（坐標變換加傳統單元細分法，其中傳統單元細分法是源點直接和單元頂點相連），得到的數值結果如表1所示。

當四邊形單元的邊長比為5時，給定4個頂點為(0.0,0.0,0.0)、(5.0,0.0,0.0)、(5.0,1.0,0.0)和(0.0,1.0,0.0)，計算結果如表2所示。

表2 四邊形單元奇異積分數值結果（邊長比5:1）

當四邊形單元的邊長比為10時，給定4個頂點為(0.0,0.0,0.0)、(10.0,0.0,0.0)、(10.0,1.0,0.0)和(0.0,1.0,0.0)，計算結果如表3所示。

從表1～表3可以看出，在高斯點數相當的前提下，隨著積分單元邊長比的增加（即單元變得狹長），本文方法的計算精度一直保持較高。

表1 四邊形單元奇異積分數值結果（邊長比1:1）

表3 四邊形單元奇異積分數值結果（邊長比10:1）

3.2 薄型平板結構

把本文方法應用到彈性薄板結構的邊界元法求解中，考慮圖5這樣一個薄板，長寬高分別為l=10 mm，l=10 mm，h=1 mm。統一量綱下彈性模量和泊松比分別為E=1 Mpa，v=0.25，質量密度為ρ=1.14。為了方便對比，在薄板邊界上施加具有解析解的二次位移場，表達式如式(13)所示。

圖5 薄板模型

把薄板模型離散成140個四邊形單元，側面用細長的四邊形單元進行離散。在x=10、z=0.5這條線上選擇一系列樣本點，采用本文方法和傳統方法的計算結果如圖6和圖7所示。

圖6 x方向的面力結果

圖7 米塞斯應力結果

由此看出，相比傳統方法（坐標變換加上傳統單元細分法，其中單元細分法是源點直接和單元頂點相連接的方法），本文方法的計算結果和解析解吻合得更好，精度更高，進一步驗證了本文方法的精確性和有效性。

4 結論

本文從彈性薄型問題邊界積分方程中存在的弱奇異積分的性質出發，提出了一種單元細分結合坐標變換的方法。該方法實現起來更加簡單，無論源點在單元的什么位置，無論單元形狀怎樣，都可以得到精確穩定的計算結果。然后通過C++編程軟件，結合UG建模和仿真平臺實現本文方法，并將其應用到薄板結構的邊界元法求解中，獲得了較好的計算結果，進一步驗證了本文方法的精確性和有效性。