無額外自由度廣義有限元的近不可壓彈-塑性分析

馬今偉, 段慶林

(1. 大連理工大學 工業裝備結構分析優化與CAE軟件全國重點實驗室, 遼寧 大連 116024;2. 大連理工大學 大連理工大學白俄羅斯國立大學聯合學院, 遼寧 大連 116024)

0 引 言

近或完全不可壓分析中的體積自鎖一直是有限元分析難題[1],它具有多種觸發形式,例如線彈性材料中無限接近0.5的Poisson比、橡膠等超彈性材料中趨于無窮的體積模量以及金屬等材料塑性應變的流動方向都會引起有限元分析中的體積自鎖．從插值的角度解釋:常規有限元的簡單形函數難以做到不違背體積近似不變約束的同時,正確反應材料和結構的變形響應．解決體積自鎖的思路大體上可以分為兩種:第一種是放松約束,最典型的代表是混合格式中的u-p方法[2];以及純位移格式中的選擇減縮積分[3]和平均體應變[4]方法．第二種思路是采用更加豐富的近似函數,如高階單元中的8節點四邊形單元和20節點六面體單元;以及p型方法中的等幾何分析、無單元Galerkin法、h-p云團法等,廣義有限元就是p型方法的一種．

廣義有限元是有限元的自然延伸,它通過引入強化函數極大地豐富了常規有限元的近似空間,并廣泛應用于裂紋擴展[5-6]、并行計算[7]等問題的模擬．傳統的廣義有限元中,強化函數的構造需要引入額外自由度,這不僅擴大了問題的求解規模,額外自由度同標準自由度之間的線性相關性還會導致剛度陣奇異．針對該問題,目前已發展了多種方法,如穩定廣義有限元[8]、正交廣義有限元[9]等．Tian[10]提出了無額外自由度的強化函數,它的構造僅需標準自由度,從而消除了線性相關性問題．這種方法與常規有限元法具有相同的計算流程,且不增加求解規模,差別僅在于近似函數不同．Xiao等[11]將這種方法推廣到了擴展有限元,并應用于裂紋擴展問題的研究．Ma等[12-13]將其推廣到了彈塑性固體的大變形分析中．在這些應用中,無額外自由度廣義有限元展現出了計算精度高、收斂穩定、網格扭曲健壯等方面的優點．

大量數值分析表明,采用高階近似函數能有效提高其近不可壓分析的能力,如:Elguedj等[14]采用不同階次的等幾何單元進行了線彈性、超彈性以及塑性材料的近不可壓分析;Chen等[15]采用再生核粒子法對橡膠梁進行了大變形分析以及塑性桿的頸縮分析等．本文的主要研究內容就是探究無額外自由度廣義有限元的近不可壓分析能力．

1 無額外自由度廣義有限元插值函數

無額外自由度廣義有限元對標量函數的插值可以寫為

(1)

其中,E為點x所在單元,NI為標準有限元形函數,uI(x)為定義在節點I處的局部近似函數,表示為

(2)

(3)

圖1 節點patch示意圖Fig. 1 Schematic diagram of patches

其中,p(x)為基底向量,它的常數形式以及二階和三階形式分別寫為

(4)

(5)

(6)

2 非線性廣義有限元節點內力線性化

控制方程弱形式的離散化可以寫為

(7)

(8)

(9)

3 近不可壓線彈性小變形分析

線彈性小變形分析忽略幾何效應,且本構模量為常量．式(9)中考慮

(10)

其中,μ和λ為Lamé常數,后者與彈性模量E及Poisson比υ的關系可寫為λ=υE/((1+υ)(1-2υ)),與體積模量κ的關系為κ=λ+2μ/3．當υ趨近于0.5時,λ和κ趨近于無窮,這意味著結構幾乎只能發生純剪切變形,而體積不能發生膨脹和壓縮．

采用經典的Cook膜問題考察無額外自由度廣義有限元在近不可壓線彈性分析中的表現．Cook膜算例的描述可參見文獻[14]．分別采用常數形式以及式(4)中二次和三次的基底向量來考察插值函數的階次對緩解體積自鎖的效果．為方便闡述,我們將采用這3種基底的方法分別表述為constant-base GFEM、quadratic-base GFEM和cubic-base GFEM．Poisson比設置為υ=0.499 9,計算網格和3種方法計算得到的變形及y方向位移云圖如圖2所示．

圖2 Cook膜算例中GFEM采用3種基向量計算的位移場和變形Fig. 2 Displacement fields and deformations obtained with the GFEM for constant, quadratic and cubic bases in the Cook membrane exmaple

采用不同密度的4種網格考察3種基底的收斂性．這些網格厚度方向始終保持為3個單元,另兩個方向的單元數分別為3×3, 6×6, 9×9, 15×15．粗略來說,右上角頂點處的位移的參考解在7.9～8.0之間,收斂性測試結果如圖3所示．

圖3 Cook膜算例中GFEM采用3種基向量在不同密度網格下收斂性測試結果Fig. 3 Convergence results of the GFEM for constant, quadratic and cubic bases in the Cook membrane example

4 近不可壓超彈性大變形分析

橡膠是典型的超彈性材料,其本構方程通過聯系第二PK應力Sij和Green應變張量Ekl給出:

(11)

(12)

其中,FiI為變形梯度．這里考察超彈性橡膠梁的純彎曲,算例的描述及Mooney-Rivlin橡膠可參見文獻[12]．采用隱式分析流程,將梁的彎矩載荷均分為100步施加,并將體積模量分別設置為κ=104,105,106,107．圖4和圖5分別展示了FEM和GFEM在第1,10,30,60,100個載荷步中梁的變形過程．

圖4 FEM在純彎曲橡膠梁中的變形Fig. 4 Deformations obtained with the FEM in the pure bending rubber beam example

圖5 GFEM在純彎曲橡膠梁中的變形Fig. 5 Deformations obtained with the GFEM in the pure bending rubber beam example

5 近不可壓J2塑性分析

J2塑性的屈服函數可寫為

(13)

(14)

(15)

其中,λ為塑性乘子,n為塑性流動方向．

塑性材料的變形過程中,塑性應變遠大于彈性應變,而塑性應變流動只能沿著偏應變方向進行,這意味著塑性變形是純剪切變形,體積無法發生膨脹或壓縮．需要說明的是,塑性變形與加載路徑相關,為保證收斂及避免可能出現的偽卸載現象需要采用正割的Newton-Raphson迭代法．此外,屈服面附近發生彈性到塑性的突變,需要采用基于應力更新算法的一致切線模量來提高收斂性．正割的Newton-Raphson迭代方法可參見文獻[16]第6章第4節,應力更新算法及一致切線模量的推導可參見文獻[17]第3章．

考察圓桿的頸縮問題,算例描述可參見文獻[15]．桿的變形集中在局部從而發生頸縮現象,而想要模擬出頸縮現象需要較為精確的內力線性化和較強的近不可壓分析能力．圖6展示了計算網格以及FEM和GFEM計算的變形及位移場,圖7展示了兩種方法得到的“頂端位移-頸縮比例”曲線和“頂端位移-支反力”曲線同試驗結果的對比,試驗數據來源于文獻[15]．

(a) 網格 (b) FEM計算的變形及位移場 (c) GFEM計算的變形及位移場 (a) The mesh (b) The FEM in the necking bar example(c) The GFEM in the necking bar example

圖7 圓桿頸縮算例中FEM和GFEM數值結果與試驗結果[15]的對比Fig. 7 Comparison between experimental data[15] and numerical results in the necking circular bar example

6 結 論

本文在彈性、超彈性及塑性不可壓分析中,得出了如下結論:

1) 彈性問題的不可壓性體現在Poisson比接近0.5時．在Cook膜算例中可以看出,采用常量形式基底向量的廣義有限元,也就是常規有限元,存在非常嚴重的體積自鎖問題,體積自鎖的緩解效果隨著插值階次的提高而變得更好．

2) 超彈性材料的不可壓性體現在超大的體積模量,常規有限元隨著體積模量的增大表現出愈發嚴重的體積自鎖,而廣義有限元方法能有效緩解這種體積自鎖．

3) 塑性材料的不可壓性體現在塑性流動只能沿著偏應變方向,導致塑性變形中體積不會膨脹和壓縮．在圓桿頸縮算例中,常規有限元無法模擬頸縮現象,而采用二次基底的廣義有限元可以．

致謝本文作者衷心感謝大連理工大學中白國際合作基金(ICR2203)對本文的資助．