郭壯志

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

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

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

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

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

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

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

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

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

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

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

2.1 求解奇異積分的單元細(xì)分技術(shù)

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

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

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

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

圖3 細(xì)長四邊形單元的細(xì)分方法

2.2 弱奇異性消除方法

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

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

圖4 三角形積分塊的坐標(biāo)變換

坐標(biāo)變換的具體過程如下：

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

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

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

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

3 數(shù)值算例

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

其中，In為數(shù)值計算結(jié)果，Ie為相應(yīng)的參考解。參考解通過使用單元細(xì)分和坐標(biāo)變換法，采用大量高斯積分點得到。為了方便計算，假設(shè)在式(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。采用本文方法和傳統(tǒng)方法（坐標(biāo)變換加傳統(tǒng)單元細(xì)分法，其中傳統(tǒng)單元細(xì)分法是源點直接和單元頂點相連），得到的數(shù)值結(jié)果如表1所示。

當(dāng)四邊形單元的邊長比為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)，計算結(jié)果如表2所示。

表2 四邊形單元奇異積分?jǐn)?shù)值結(jié)果（邊長比5:1）

當(dāng)四邊形單元的邊長比為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)，計算結(jié)果如表3所示。

從表1～表3可以看出，在高斯點數(shù)相當(dāng)?shù)那疤嵯拢S著積分單元邊長比的增加（即單元變得狹長），本文方法的計算精度一直保持較高。

表1 四邊形單元奇異積分?jǐn)?shù)值結(jié)果（邊長比1:1）

表3 四邊形單元奇異積分?jǐn)?shù)值結(jié)果（邊長比10:1）

3.2 薄型平板結(jié)構(gòu)

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

圖5 薄板模型

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

圖6 x方向的面力結(jié)果

圖7 米塞斯應(yīng)力結(jié)果

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

4 結(jié)論

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