潘先云，周楓林，王瑾元，袁小涵，欽 宇

（湖南工業(yè)大學(xué) 機(jī)械工程學(xué)院，湖南 株洲 412007）

1 研究背景

結(jié)構(gòu)靜彈性問題研究結(jié)構(gòu)在受到外力作用下而產(chǎn)生的應(yīng)力、應(yīng)變和位移的變化，以維持結(jié)構(gòu)在合理變化范圍之內(nèi)，保證應(yīng)用場(chǎng)景的安全性，為工程設(shè)計(jì)提供參考。如裴智毅等[1]運(yùn)用彈性力學(xué)對(duì)鋼橋面鋪裝層進(jìn)行了應(yīng)力分析，為鋼橋面STC（super-through concrete）鋪裝的設(shè)計(jì)和箱梁橫隔板截面尺寸的選取提供了參考；郭一竹等[2]研究了彈性伸桿的力學(xué)性能，并且對(duì)其剛度進(jìn)行了分析。結(jié)構(gòu)靜彈性問題分析的數(shù)值方法以有限元法（finite element method，F(xiàn)EM）[3-5]為代表，其中心思想是將連續(xù)的彈性體離散為有限個(gè)單元，再根據(jù)彈性力學(xué)基本方程和變分原理建立單元節(jié)點(diǎn)間位移與力的關(guān)系，形成單元?jiǎng)偠染仃嚕缓髮⑺袉卧獎(jiǎng)偠染仃囌沓煽傮w剛度矩陣，最后通過引入邊界條件，生成關(guān)于未知節(jié)點(diǎn)位移的線性方程組。ABAQUS、ANSYS等有限元法仿真軟件的出現(xiàn)，使有限元法變得簡(jiǎn)單高效，并在機(jī)械結(jié)構(gòu)分析、土木工程領(lǐng)域得到了廣泛應(yīng)用。但是，有限元法仍存在不足之處。首先，該方法需要將整個(gè)結(jié)構(gòu)區(qū)域進(jìn)行離散，處理的數(shù)據(jù)量較大，計(jì)算時(shí)間較長(zhǎng)；其次，難以精確處理無限域和應(yīng)力集中等問題，上述缺陷影響了有限元法的精度和實(shí)際應(yīng)用范圍。因此，研究者在有限元法的基礎(chǔ)上發(fā)展了一種新的數(shù)值計(jì)算方法：邊界元法[6-8]（boundary element method，BEM）。相比于有限元法，邊界元法只需要對(duì)邊界進(jìn)行離散，不涉及域內(nèi)，因而將問題的維數(shù)降低，擁有較高的計(jì)算效率和計(jì)算精度。邊界元法在彈性力學(xué)中有許多研究成果和成功案例。如曾華等[9]利用雙互易邊界元法求解了彈性問題，證明其具有較高的精度；徐俊祥等[10]采用直接邊界元法對(duì)彈性體整個(gè)域內(nèi)的位移場(chǎng)和應(yīng)力場(chǎng)進(jìn)行了求解。邊界元法在其它領(lǐng)域也有著廣泛應(yīng)用。

本文運(yùn)用雙互易邊界元法（dual reciprocity boundary element method，DRBEM）分析結(jié)構(gòu)靜彈性中含有體力項(xiàng)的非齊次問題，以及在無解析解的實(shí)際工況下，通過雙互易邊界元法和有限元法分析結(jié)果進(jìn)行對(duì)比，驗(yàn)證雙互易邊界元法在分析結(jié)構(gòu)靜彈性問題上的有效性和精確性，并為處理非齊次彈性力學(xué)問題提供一種新的方法。有學(xué)者曾對(duì)邊界元法和有限元法進(jìn)行了對(duì)比研究，比如豆中強(qiáng)[11]利用有限元分析軟件——ANSYS軟件和邊界元法求解二維彈性靜力學(xué)問題，并將求解得到的結(jié)果與精確解進(jìn)行對(duì)比；周楓林等[12]分別用有限元法和邊界元法分析了印刷機(jī)滾筒結(jié)構(gòu)剛度，證明了邊界元法是用于結(jié)構(gòu)剛度分析的一種高效方法。與上述研究的不同之處在于，本研究采用指數(shù)型基函數(shù)作為雙互易邊界元法的插值函數(shù)，對(duì)結(jié)構(gòu)進(jìn)行靜彈性力學(xué)分析，驗(yàn)證該方法的有效性，并將DRBEM和FEM對(duì)比，驗(yàn)證DRBEM在實(shí)際工況下結(jié)構(gòu)靜彈性分析中的精確性。

2 靜彈性力學(xué)雙互易邊界元法

2.1 靜彈性力學(xué)的理論基礎(chǔ)

課題組采用雙互易邊界元法，分析正方體模型在受外力載荷條件下引起的彈性變形情況。

三維靜彈性力學(xué)問題的基本方程如下：

平衡微分方程為

幾何方程為

物理方程（本構(gòu)方程）為

式（1）～（3）中：i,j,k=1, 2, 3；

σ為應(yīng)力；

b為體力；

u為位移；

δ為Kronecker函數(shù)；

G為剪切模量；

λ為拉梅常數(shù)，且

其中，E為彈性模量，v為泊松比；

ε為應(yīng)變，且

將式（2）代入式（3）中，得到如下應(yīng)力分量用位移分量的表示：

再將式（6）代入式（1）中，可得到如下用位移表示的平衡方程（Navier方程）：

2.2 靜彈性力學(xué)雙互易邊界元法描述

上文得到用位移表示的平衡方程如下：

位移邊界條件為

面力邊界條件為

式中：n為邊界外法向量；

Ω為彈性體所占區(qū)域；

Г1和Г2分別為位移邊界和面力邊界；

可用加權(quán)余量法、功的互等定理等建立邊界積分方程，現(xiàn)利用加權(quán)余量法、散度定理將式（1）、（9）和（10）建立彈性力學(xué)邊界積分方程，簡(jiǎn)化后如下：

其中，r為場(chǎng)點(diǎn)與源點(diǎn)之間的距離，l=1, 2, 3；

式中，θ為源點(diǎn)處的立體角，

式（11）是在bj=0情況下的邊界積分方程，也就是計(jì)算域內(nèi)沒有外力作用，式（12）存在體力項(xiàng)。式（11）和（12）描述的是已知邊界Г上位移uk和面力pk計(jì)算區(qū)域內(nèi)所有點(diǎn)p在k方向上的位移。

對(duì)于邊界積分方程式（12）中含域內(nèi)體力積分非齊次項(xiàng)問題，利用雙互易方法處理。用徑向基函數(shù)對(duì)式（12）中的體力項(xiàng)進(jìn)行插值擬合，插值表達(dá)式為

式中：φj為插值函數(shù)；

H和L分別為邊界點(diǎn)和內(nèi)部點(diǎn)的數(shù)量。

值得注意的是，若體力項(xiàng)函數(shù)未知，則徑向基函數(shù)插值系數(shù)也將是未知的。找到徑向基函數(shù)對(duì)應(yīng)式（7）特定的位移u滿足：

利用分部積分可以將式（12）轉(zhuǎn)換為完全邊界積分方程：

m=1, 2, 3。

邊界積分方程式（11）和（18）只有邊界積分，所以可以通過離散邊界進(jìn)行求解計(jì)算，將位移和面力離散為如下形式：

M為每一個(gè)單元的節(jié)點(diǎn)數(shù)量。

將計(jì)算邊界Г離散為N個(gè)邊界單元，再將式（19）和（20）代入式（11）和（18）中，得到離散形式的邊界積分表達(dá)式如下：

式中：Γf為第f個(gè)單元的邊界。

最終離散后式（21）和式（22）可形成如下線性代數(shù)方程組：

將所有未知量移到方程的左邊并整理可得

通過式（25）所示代數(shù)方程組可以求出所有邊界節(jié)點(diǎn)的未知量，包括面力和位移。

2.3 指數(shù)型插值函數(shù)及其特解

本文采用指數(shù)型插值函數(shù)作為雙互易邊界元法的徑向基函數(shù)，對(duì)體力項(xiàng)進(jìn)行插值擬合，指數(shù)型插值函數(shù)表達(dá)式如下：

式中：c為形狀參數(shù)；

r為場(chǎng)點(diǎn)Q到源點(diǎn)P之間的距離，且

對(duì)應(yīng)的指數(shù)型插值函數(shù)在彈性問題上的位移特解表達(dá)式為

式中：

其中：

上述位移特解并不是唯一的。該位移特解是通過Papkovich勢(shì)函數(shù)法推導(dǎo)出來的滿足平衡方程的一種形式，并且保證了g(r)與q(r)在全域內(nèi)連續(xù)且有一階導(dǎo)數(shù)。不同的RBF對(duì)應(yīng)的特解也會(huì)改變，將位移特解導(dǎo)入胡克定律與平衡方程，同樣可以求解出應(yīng)變、應(yīng)力和位移特解。

3 數(shù)值算例

為了驗(yàn)證雙互易邊界元法在靜力學(xué)分析上的有效性和準(zhǔn)確性，設(shè)計(jì)了兩個(gè)算例。第一個(gè)算例是有解析解情況下，通過DRBEM方法得出數(shù)值解與解析解對(duì)比，驗(yàn)證雙互易邊界元法的有效性。第二個(gè)算例是在無解析解的實(shí)際工況下，將DRBEM結(jié)果與有限元分析軟件結(jié)果進(jìn)行對(duì)比，進(jìn)一步驗(yàn)證雙互易邊界元法的有效性和準(zhǔn)確性。

根據(jù)對(duì)指數(shù)型基函數(shù)作為雙互易邊界元法的徑向基函數(shù)的分析研究，當(dāng)指數(shù)型基函數(shù)的形狀參數(shù)c=0.2時(shí)，數(shù)值計(jì)算精度最高，所以本文算例雙互易邊界元法中，指數(shù)型徑向基函數(shù)形狀參數(shù)c取0.2。

3.1 球體邊界元靜力學(xué)分析

本算例研究半徑為2 mm的球體，笛卡爾坐標(biāo)中心為球心。球體的材料參數(shù)設(shè)置如下：密度為1.14 t/mm3；泊松比為0.25；楊氏模量為1 000 MPa。通過雙互易邊界元法對(duì)球體進(jìn)行彈性力學(xué)分析，其網(wǎng)格劃分如圖1所示，網(wǎng)格數(shù)量為2 250個(gè)。

圖1 球體網(wǎng)格劃分Fig.1 Sphere meshing

通過改變網(wǎng)格大小、增加網(wǎng)格數(shù)量，計(jì)算得到節(jié)點(diǎn)上面力以及域內(nèi)應(yīng)力的數(shù)值解，將其與精確解對(duì)比得到相對(duì)誤差，精確解是通過控制方程式（8）得到的，其中相對(duì)誤差定義為

式中：v為相關(guān)變量；

為精確解；

為數(shù)值解；

|v|max為精確解中的最大值；

N為節(jié)點(diǎn)數(shù)量。

在邊界元法中，3個(gè)位移邊界條件如下：

將式（29）代入式（6）中可以得到應(yīng)力的解析解：

圖2和圖3所示結(jié)果是利用雙互易邊界元法對(duì)球體進(jìn)行靜彈性力學(xué)分析得到的，其中圖2所示為單元數(shù)量對(duì)面力相對(duì)誤差的影響，圖3所示為單元數(shù)量對(duì)域內(nèi)應(yīng)力的影響。

圖2 單元數(shù)量對(duì)面力相對(duì)誤差的影響曲線Fig.2 Influence curves of element number on relative error of surface force

圖3 單元數(shù)量對(duì)域內(nèi)應(yīng)力相對(duì)誤差的影響曲線Fig.3 Influence curves of element number on relative error of stress in domain

從圖2和圖3中可以看到，隨著單元數(shù)量增多，節(jié)點(diǎn)數(shù)量增加，節(jié)點(diǎn)上各面力相對(duì)誤差值和域內(nèi)應(yīng)力的相對(duì)誤差值均隨之減小，數(shù)值解均收斂于解析解，從而得到穩(wěn)定面力和應(yīng)力結(jié)果。這一結(jié)果說明，采用雙互易邊界元法分析球體的靜彈性力學(xué)問題的方法是有效的。

3.2 DRBEM與FEM對(duì)比分析

本算例在實(shí)際工況下，利用DRBEM和FEM方法，對(duì)立方體模型進(jìn)行靜彈性力學(xué)分析，以驗(yàn)證實(shí)際工況下DRBEM分析結(jié)構(gòu)靜彈性力學(xué)方面的準(zhǔn)確性和可行性。所用模型為邊長(zhǎng)為20 mm的立方體，笛卡爾坐標(biāo)中心為立方體的質(zhì)心位置，其坐標(biāo)如表1所示。

表1 立方體取樣點(diǎn)坐標(biāo)Table 1 Cube sampling point coordinates

取樣點(diǎn)位置如圖4所示。立方體材料參數(shù)設(shè)置如下：密度0.007 9 t/mm3；泊松比為0.28；楊氏模量為2 000 MPa。邊界條件設(shè)置如下：左表面施加壓強(qiáng)為100 MPa；給定Z方向體力為-20 t/(s2·mm2)；下表面固定約束。

圖4 正方體取樣點(diǎn)位置分布圖Fig.4 Location distribution of cube sampling points

邊界元網(wǎng)格數(shù)量為1 608，采用非連續(xù)三角形單元類型，網(wǎng)格模型見圖5。有限元網(wǎng)格數(shù)量為24 389，采用線性六面體單元類型，網(wǎng)格模型如圖6所示。利用有限元方法仿真軟件對(duì)立方體進(jìn)行靜力學(xué)分析，得到的位移云圖如圖7所示。

圖5 立方體邊界元網(wǎng)格模型Fig.5 Cube boundary element mesh model

圖6 立方體有限元網(wǎng)格模型Fig.6 Cube finite element mesh model

圖7 有限元位移云圖Fig.7 Finite element displacement nephogram of cube model

圖8所示為以邊界元法和有限元法的位移結(jié)果比較圖，表2給出了求得的取樣點(diǎn)的位移結(jié)果具體數(shù)據(jù)。

圖8 DRBEM和FEM位移結(jié)果比較Fig.8 Comparison of displacement results between DRBEM and FEM

綜合圖8和表2中DRBEM和FEM的取樣點(diǎn)位移數(shù)值情況看，兩者計(jì)算位移數(shù)值結(jié)果高度擬合，該算例證明了在實(shí)際工況下雙互易邊界元法在分析結(jié)構(gòu)靜彈性力學(xué)方面的準(zhǔn)確性和可行性。

表2 DRBEM和FEM的取樣點(diǎn)位移結(jié)果Table 2 Displacement results of DRBEM and FEM sampling points mm

4 結(jié)語

本文介紹了雙互易邊界元法在靜彈性力學(xué)問題中的應(yīng)用，以及雙互易邊界元法與有限元法在靜彈性力學(xué)分析中的對(duì)比。算例結(jié)果表明，雙互易邊界元法是分析結(jié)構(gòu)靜彈性力學(xué)問題的一種有效方法。在實(shí)際工況下，通過DRBEM與FEM在靜彈性力學(xué)分析上的對(duì)比，進(jìn)一步證明DRBEM方法在分析結(jié)構(gòu)靜彈性力學(xué)方面的準(zhǔn)確性和可行性，并且具有精度高等特點(diǎn)。由此可見，DRBEM為分析結(jié)構(gòu)靜彈性力學(xué)問題提供了一種新的方法，并可以將雙互易邊界元法應(yīng)用到解決其他領(lǐng)域含域積分項(xiàng)的非齊次問題。