馬今偉，段慶林，陳嵩濤

(大連理工大學 工業(yè)裝備結構分析國家重點實驗室，大連 116024)

1 引 言

廣義有限元法[1]能夠方便地引入待求解問題相關的局部強化函數，具有計算精度高、易于程序實現(對已有的有限元程序改動較小)等優(yōu)點，引發(fā)了廣泛關注。該方法一般認為源于單位分解法[2]和單位分解有限元法[3]，甚至還可追溯到更早期的研究[4,5]。由于局部強化函數可以根據具體問題方便靈活地引入，廣義有限元法已廣泛應用于各類問題，如裂紋擴展[6]、形狀優(yōu)化[7]、非光滑界面[8]和水力壓裂[9]等。

在廣義有限元法[1]中，局部強化函數的引入依賴于額外自由度。以本文考慮的靜力問題為例，每個計算節(jié)點除標準的位移自由度之外，還需設置與局部強化函數數目相等的額外自由度。這導致了求解規(guī)模的急劇增大。更為嚴重的是，依賴額外自由度引入局部強化函數會造成線性相關性問題[3,10,11]，導致剛度陣奇異。針對該問題，目前已有一些行之有效的方法，如穩(wěn)定廣義有限元方法(SGFEM)[12]以及正交廣義有限元方法(OGFEM)[13]等。這些方法能有效消除線性相關性問題，但仍然需要使用額外自由度，不能減小求解規(guī)模，實現起來也與標準有限元法有較大差異。

田榮[14]采取了一種不同的思路，提出不使用額外自由度引入局部強化函數，建立了無額外自由度廣義有限元法。該方法不僅消除了前述線性相關問題，而且節(jié)點自由度與標準有限元法相同，同時仍然保留了廣義有限元方法的諸多優(yōu)勢，如局部強化函數仍然能夠方便靈活地引入?；谶@一思想，田榮等[15,16]也對廣泛應用于裂紋問題的擴展有限元法進行了改進，建立了iXFEM(improved XFEM)方法。同樣地，無額外自由度廣義有限元法作為對原GFEM方法的重要改進，在本文中簡記為iGFEM。

目前，對于本文考慮的靜力問題，iGFEM仍然僅限于線彈性小變形分析。本文將基于考慮幾何和物理非線性的計算方法理論[17-20]，將該方法推廣至彈塑性大變形問題的非線性分析。雖然，有限元法的非線性分析相對成熟，并已發(fā)展了諸多計算軟件，包括我國自主軟件SiPESC[21]。但是，低階單元非線性分析的計算精度往往難以令人滿意。iGFEM 方法可使用低階單元的計算網格建立高階的插值近似，若能有效改善非線性分析的計算精度，將具有重要的研究意義和潛在的應用價值。

2 控制方程

參考初始構型，靜力平衡方程及邊界條件為

?Pi j/?Xj+bi=0inΩ0

(1)

(2)

Wext-Wint=0

(3)

式中Wint和Wext分別表示內力和外力虛功:

(4)

(5)

式中δui為虛位移。弱形式(3)的進一步空間離散依賴于位移場的插值近似，這是iGFEM方法的核心內容。

3 無額外自由度的廣義有限元近似

廣義有限元方法的位移近似可寫為

(6)

(7)

該式可進一步寫為

(8)

式中

(9)

圖1 計算網格上的節(jié)點片

(10)

其中

(11)

式中ω(x)為MLS的權函數，pK=p(xK)，本文取為如下的二次基底

(12)

4 節(jié)點力及其線性化

為得到離散化方程，將iGFEM的位移近似式(8)代入弱形式(3)，并通過位移的獨立變分可得到

(13)

(14)

(15)

式(13)是以位移為基本未知量的非線性方程，需采用Newton-Raphson迭代求解。為此，對節(jié)點內力作如下線性化

(16)

式(16)的推導利用了以下關系：

(17)

(18)

其中

(19)

(20)

應說明的是，本文僅考慮不隨構形變化的外載，因而無須對節(jié)點外力進行線性化。

5 本構模型

本文考慮的亞彈-塑性材料的彈性本構關系可寫為

(21)

(22)

(23)

(24)

式中γ為塑性參數，rk l為塑性流動方向。對于本文考慮的J2關聯塑性有rk l=?f/?σk l，f為屈服面。

(25)

式中q為內變量，在本文僅考慮為等效塑性應變，‖ξ‖為等效應力，K(α)=σY+Hα為后繼屈服強度，σY為初始屈服強度，α為等效塑性應變，H為塑性模量。將式(23，24)代入式(21)可得

(26)

(γ≥0，f≤0)(27)

對于本文考慮的超彈性材料，本構模型可由第二PK應力Si j和格林應變Ek l寫為

(28)

(29)

將上述兩種本構模型代入式(18)，則節(jié)點內力率可寫為

(30)

(31)

(32)

6 數值算例

采用三個算例考察本文發(fā)展的非線性iGFEM 方法處理彈塑性大變形問題的有效性和計算精度。為作比較，對標準線性有限元方法也進行了數值實現，并標記為FEM。

6.1 淺拱

如圖2所示，兩端固支弧形淺拱在頂部中心點A處受向下的集中力作用。淺拱中性軸半徑為100 mm，拱厚為2 mm，淺拱對應圓心角為28.06°，材料為超彈性，楊氏模量E=4.8×103N/mm2，泊松比υ=0。iGFEM方法計算得到的淺拱變形及相應的σx x應力場如圖3所示。圖4比較了兩種方法的載荷-位移曲線。在相同的計算網格下，iGFEM 的計算結果與解析解吻合，FEM則發(fā)生了明顯的偏離。應說明的是，本文FEM僅指標準的線性三角形單元，若采用高階單元或使用精巧的單元技術，有限元法的計算精度將會得到顯著提高。

圖2 受集中力作用的淺拱

圖3 淺拱變形及σxx應力場

圖4 淺拱A點處的載荷-位移曲線

6.2 懸臂梁

如圖5所示，懸臂梁一端固支，另一端受向上大小為F=10 N的集中力作用，梁長為20 mm，厚度為1 mm，材料為超彈性，楊氏模量E=4.8×103N/mm2，泊松比υ=0。iGFEM計算得到的懸臂梁變形過程如圖6所示。圖7比較了懸臂梁自由端的載荷-撓度曲線，再次展現了iGFEM方法的高精度。

6.3 圓桿頸縮

采用圓桿頸縮問題考察本文方法對于亞彈-塑性材料大變形計算的有效性。如圖8所示，圓桿長為53.334 mm，半徑為6.413 mm，兩端各施加 7 mm 的固定位移。由于對稱性，只取圓桿的1/4

圖5 自由端受集中力的懸臂梁

圖6 懸臂梁變形過程

圖7 懸臂梁載荷-撓度響應對比

進行數值模擬。在圓桿中間截面處引入幾何缺陷，該處的半徑為兩端半徑的98.2%。具體的材料參數和實驗數據可以參考文獻[22]。

采用圖9所示的兩種計算網格，FEM和iGFEM 計算得到的頸縮曲線與實驗數據[22]在 圖10 進行了比較。兩種方法計算得到的圓桿頸縮變形分別如圖11和圖12所示?？梢钥闯觯現EM粗細兩種網格的計算結果有明顯差異，而iGFEM粗細網格的計算結果基本保持一致，這在一定程度上驗證了iGFEM方法在粗網格下具有良好的高精度特性。

圖8 圓桿頸縮算例

圖9 圓桿頸縮算例的計算網格

圖10 圓桿頸縮算例的頸縮-位移曲線

圖11 FEM方法得到的圓桿頸縮變形

圖12 i GFEM方法得到的圓桿頸縮變形

7 結論與展望

本文建立了幾何和物理非線性分析的無額外自由度廣義有限元方法。數值結果表明，本文方法不僅能有效分析彈塑性大變形問題，而且展現出優(yōu)于傳統(tǒng)線性有限元方法的計算精度。后續(xù)工作將把該方法推廣至三維，面向復雜工程問題，有力推動彈塑性大變形高精度計算方法的發(fā)展。