曾 華， 周楓林， 余江鴻

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

0 引言

在工程實(shí)踐中往往會遇到一些彈性力學(xué)方面的問題，而一般的彈性力學(xué)問題用解析方法求解就會非常麻煩，此時(shí)，采用數(shù)值方法來解決這類問題更為有效[1]。 到目前為止，人們所熟悉的工程數(shù)值方法主要有：有限單元法[2]、Ritz 法[3]、Galerkin 法[4]以及邊界元法[5]。 近年來，隨著計(jì)算機(jī)技術(shù)飛速發(fā)展，邊界元法得到了廣泛應(yīng)用，邊界元法最重要的特點(diǎn)是降低維度，從而提高計(jì)算效率。由于邊界元法具有較高的計(jì)算精度， 在處理無限域問題或應(yīng)力集中問題的精度比有限元分析要好， 所以這種方法受到了國內(nèi)外許多研究者的重視。 邊界元法作為一種常用的數(shù)值分析方法，在彈性力學(xué)、流體力學(xué)、熱傳導(dǎo)、聲場，以及電磁場等研究領(lǐng)域得到了廣泛應(yīng)用[6]。 在邊界元法中，由于邊界積分方程對于含體力項(xiàng)的彈性問題求解起來十分困難，需要采用一種更有效的方法來解決體力項(xiàng)的問題，而雙重互易法（DRM）的出現(xiàn)就很好地解決體力項(xiàng)的問題，近年來雙重互易法被許多學(xué)者所研究，其中文獻(xiàn)[7-10]研究較多， 雙重互易法實(shí)現(xiàn)了將體力項(xiàng)產(chǎn)生的域內(nèi)積分轉(zhuǎn)化為等效的邊界積分， 從而有效避免了域內(nèi)積分的復(fù)雜計(jì)算， 本文通過編寫邊界元程序驗(yàn)證了雙重互易法在求解彈性問題的有效性。

1 彈性力學(xué)問題的邊界元法

描述線彈性問題的平衡方程、幾何方程和本構(gòu)方程為：

式中：i，j—方向；σij—應(yīng)力；εkl—應(yīng)變；bi—體積力；Eijkl—彈性模量；uij—位移，uij，j表示uij在j 方向上的導(dǎo)數(shù)；V—域內(nèi)體積。 若要求出方程的解還需要包含一定的邊界條件，可以寫成：

式中：pi—面力；nj—邊界外法線方向的法向量；S—邊界面積。

通過Betti 功互等定理（簡稱Betti 定理）我們可以建立邊界積分方程，形如：

式中：c—光滑系數(shù)，一般光滑邊界點(diǎn)c=1/2， u*和p*分別表示位移基本解和面力基本解， 其值與節(jié)點(diǎn)位置、 泊松比、彈性模量有關(guān)。

2 雙重互易法

邊界積分方程中如果在無體力項(xiàng)的情況下， 只要邊界條件中位移或者面力中的一項(xiàng)是已知的， 就可以通過邊界積分方程求出另外一項(xiàng)， 但是由于在實(shí)際工程應(yīng)用中，物體往往都有體力項(xiàng)存在，所以對于含體力項(xiàng)的彈性力學(xué)方程求解問題，原本的邊界元方法就不太適用。

為了解決體力項(xiàng)的問題，采用了雙重互易法。即當(dāng)定義中考慮體力項(xiàng)的區(qū)域積分時(shí)， 可以將體力項(xiàng)積分近似的表示為：

式中：H 和G 矩陣分別由面力和位移基本解組成；F-1表示體積力的逆矩陣。

3 數(shù)值算例

為了驗(yàn)證雙互易邊界元法的有效性，在Visual studio的軟件環(huán)境下設(shè)計(jì)兩個(gè)邊界元法分析程序。 第一個(gè)是在給定邊界位移的條件下計(jì)算立方體模型域內(nèi)的位移變化， 第二個(gè)是計(jì)算球體模型的域內(nèi)位移變化。

3.1 立方體模型域內(nèi)位移變化

建立一個(gè)邊長為2 的立方體簡易模型， 網(wǎng)格劃分如圖1 所示，該立方體的網(wǎng)格大小為0.1，采用的邊界節(jié)點(diǎn)數(shù)目為792， 所構(gòu)建的二次三節(jié)點(diǎn)單元為1332 個(gè)， 內(nèi)部RBF 插值點(diǎn)個(gè)數(shù)為1895， 內(nèi)部插值點(diǎn)位置分布如圖2 所示，材料參數(shù)：彈性模量設(shè)置為1，泊松比0.25，材料密度1.14。 立方體的空間位置范圍為：

圖1 立方體網(wǎng)格劃分

圖2 插值點(diǎn)位置分布

其中網(wǎng)格大小與節(jié)點(diǎn)數(shù)量的關(guān)系如表2 所示，m 表示網(wǎng)格大小，N 表示節(jié)點(diǎn)數(shù)量，L 表示域內(nèi)插值點(diǎn)數(shù)量。選取5 個(gè)域內(nèi)節(jié)點(diǎn)作為參考點(diǎn)，見表1，x 方向上的計(jì)算結(jié)果如表2 所示。

表1 參考點(diǎn)坐標(biāo)

表2 網(wǎng)格大小對應(yīng)的節(jié)點(diǎn)數(shù)量

從表3 中可以看出，當(dāng)固定y，z 方向的坐標(biāo)值，隨機(jī)選擇域內(nèi)五組參考點(diǎn)， 各點(diǎn)在x 方向上的位移計(jì)算結(jié)果與解析解的吻合度比較好，說明了采用此方法是有效的，為了進(jìn)一步證明方法的準(zhǔn)確性， 驗(yàn)證了在不同網(wǎng)格大小下計(jì)算結(jié)果與解析解的變化，計(jì)算A1 在x 方向上位移相對差值如圖3 所示。

從圖3 可以看出，隨著插值點(diǎn)的不斷增加，位移相對誤差也會越來越小， 計(jì)算結(jié)果隨著網(wǎng)格大小的增加呈單調(diào)遞減的趨勢，當(dāng)網(wǎng)格大小為0.1 時(shí)，A1 的位移相對誤差為3.1393%，由此可以看出，只要選擇合適的網(wǎng)格大小，雙互易邊界元法具有良好的的精度。

表3 x 方向上的位移結(jié)果

圖3 A1 在x 分向上的位移相對差值

3.2 球體模型域內(nèi)位移變化

建立一個(gè)半徑為2 的球體模型， 網(wǎng)格劃分如圖4 所示， 該球體的網(wǎng)格大小為0.3 采用的邊界節(jié)點(diǎn)數(shù)目為581，所構(gòu)建的二次三節(jié)點(diǎn)單元為1066 個(gè)，內(nèi)部RBF 插值點(diǎn)個(gè)數(shù)為2003，內(nèi)部插值點(diǎn)位置分布如圖5 所示，材料參數(shù)與正方體相同，球體坐標(biāo)原點(diǎn)為（0，0，0）。 空間坐標(biāo)為

圖4 球體網(wǎng)格劃分

圖5 插值點(diǎn)位置分布

邊界條件：

同樣選取5 個(gè)參考點(diǎn)， 見表4，y 方向的結(jié)果如表5所示，通過表5 中計(jì)算結(jié)果與解析解進(jìn)行對比可知，此方法得到的結(jié)果與精確解吻合度很好， 兩組實(shí)驗(yàn)都反映了采用雙互易邊界元法在求解彈性問題上是有效的。

表4 參考點(diǎn)坐標(biāo)

表5 y 方向的位移結(jié)果

4 結(jié)論

運(yùn)用雙互易邊界元法避免了復(fù)雜的域內(nèi)積分計(jì)算；通過兩組數(shù)值分析， 驗(yàn)證了雙互易邊界元法在求解彈性問題上的有效性，并隨著網(wǎng)格劃分得越來越細(xì)，計(jì)算結(jié)果具有更好的精度； 不足之處在于算法并未采用任何加速計(jì)算，因此提高計(jì)算效率是進(jìn)一步需要完成的工作。