含非均勻體力機械結構彈性力學問題的雙互易邊界元法

曾 華， 周楓林， 余江鴻

（湖南工業大學 機械工程學院， 湖南 株洲412007）

0 引言

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

1 彈性力學問題的邊界元法

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

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

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

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

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

2 雙重互易法

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

為了解決體力項的問題，采用了雙重互易法。即當定義中考慮體力項的區域積分時， 可以將體力項積分近似的表示為：

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

3 數值算例

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

3.1 立方體模型域內位移變化

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

圖1 立方體網格劃分

圖2 插值點位置分布

其中網格大小與節點數量的關系如表2 所示，m 表示網格大小，N 表示節點數量，L 表示域內插值點數量。選取5 個域內節點作為參考點，見表1，x 方向上的計算結果如表2 所示。

表1 參考點坐標

表2 網格大小對應的節點數量

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

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

表3 x 方向上的位移結果

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

3.2 球體模型域內位移變化

建立一個半徑為2 的球體模型， 網格劃分如圖4 所示， 該球體的網格大小為0.3 采用的邊界節點數目為581，所構建的二次三節點單元為1066 個，內部RBF 插值點個數為2003，內部插值點位置分布如圖5 所示，材料參數與正方體相同，球體坐標原點為（0，0，0）。 空間坐標為

圖4 球體網格劃分

圖5 插值點位置分布

邊界條件：

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

表4 參考點坐標

表5 y 方向的位移結果

4 結論

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