精確計(jì)算彈性薄型結(jié)構(gòu)中奇異積分方法

郭壯志

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

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

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

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

1 薄型問(wèn)題的奇異邊界積分方程

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

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

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

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

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

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

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

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

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

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

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

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

2.2 弱奇異性消除方法

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

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

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

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

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

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

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

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

3 數(shù)值算例

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

其中，In為數(shù)值計(jì)算結(jié)果，Ie為相應(yīng)的參考解。參考解通過(guò)使用單元細(xì)分和坐標(biāo)變換法，采用大量高斯積分點(diǎn)得到。為了方便計(jì)算，假設(shè)在式(4)中，f(P,Q)=1.0。下述第一個(gè)例子考慮式（4）形式具有弱奇異性的邊界積分。

3.1 奇異積分算例

本例以四邊形單元為例來(lái)驗(yàn)證本文方法。四邊形單元的四個(gè)頂點(diǎn)為(0.0,0.0,0.0)、(1.0,0.0,0.0)、(1.0,1.0,0.0)和(0.0,1.0,0.0)，邊長(zhǎng)比為a=1。采用本文方法和傳統(tǒng)方法（坐標(biāo)變換加傳統(tǒng)單元細(xì)分法，其中傳統(tǒng)單元細(xì)分法是源點(diǎn)直接和單元頂點(diǎn)相連），得到的數(shù)值結(jié)果如表1所示。

當(dāng)四邊形單元的邊長(zhǎng)比為5時(shí)，給定4個(gè)頂點(diǎn)為(0.0,0.0,0.0)、(5.0,0.0,0.0)、(5.0,1.0,0.0)和(0.0,1.0,0.0)，計(jì)算結(jié)果如表2所示。

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

當(dāng)四邊形單元的邊長(zhǎng)比為10時(shí)，給定4個(gè)頂點(diǎn)為(0.0,0.0,0.0)、(10.0,0.0,0.0)、(10.0,1.0,0.0)和(0.0,1.0,0.0)，計(jì)算結(jié)果如表3所示。

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

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

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

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

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

圖5 薄板模型

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

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

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

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

4 結(jié)論

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