一種新型基于共用系數(shù)的二維ANCF 高階梁單元

周川，趙春花，獨(dú)亞平，郭嘉輝

（201620 上海市 上海工程技術(shù)大學(xué) 機(jī)械與汽車工程學(xué)院）

0 引言

絕對節(jié)點(diǎn)坐標(biāo)法（簡記為ANCF）[1]是由美國伊利諾伊大學(xué)的夏巴那教授提出的一種解決多體系統(tǒng)大變形大轉(zhuǎn)動的新型有限元方法。不同于牛頓歐拉法[2]、凱恩法[3]等其它的多體系統(tǒng)建模方法[4-6]，ANCF 法拋棄了傳統(tǒng)有限元小變形假設(shè)，將轉(zhuǎn)角以梯度形式表示引入到節(jié)點(diǎn)自由度中，得到的質(zhì)量矩陣為常數(shù)陣且不存在科氏力、離心力，提出至今受到越來越多學(xué)者們的關(guān)注和研究。

研究前期，Omar 等[7]提出的二維ANCF 全參數(shù)剪切梁單元通過一般連續(xù)介質(zhì)力學(xué)方法構(gòu)建單元彈性力，但一般連續(xù)介質(zhì)力學(xué)方法構(gòu)建的單元存在一些閉鎖問題[8]；沈振興[9]和張大羽等[8]學(xué)者認(rèn)為，全參數(shù)梁的位移模式在橫向橫縱向的插值不一致，導(dǎo)致無法滿足廣義胡克定律，產(chǎn)生了偽應(yīng)力，使得單元計算性能下降。為解決該問題，Zhao 等[10]通過共用系數(shù)法在橫向插值中引入y2，使得全參數(shù)剪切梁橫縱向插值保證了2 階次，取得了不錯效果；Patel 等[11]基于獨(dú)立系數(shù)引入曲率坐標(biāo)，構(gòu)建新的橫向高階單元，橫縱向插值達(dá)到2 階一致，有效減輕了閉鎖問題，同時還提出了應(yīng)變分裂法解決ANCF 梁和板單元中存在的閉鎖問題；Li 等[12]在二維全參數(shù)梁基礎(chǔ)上增加曲率自由度rxy，同時通過共用插值項系數(shù)引入了橫向高階x2y2，但是計算結(jié)果并沒有很好地改善；於祖慶[13]提出了幾種絕對節(jié)點(diǎn)坐標(biāo)單元，并將新單元用于車輛鋼板彈簧和輪胎的動力學(xué)建模中；徐齊平[14]提出了超彈性ANCF 單元，并將其用于建立氣動軟體機(jī)器人的精確動力學(xué)模型。

本文在Li 提出的單元自由度[12]基礎(chǔ)上，通過共用系數(shù)引入橫向高階插值項，提出了一種新的梁單元。通過理論推導(dǎo)結(jié)合數(shù)值實例，證明了新單元的可行性，能夠緩解泊松閉鎖問題，提高單元收斂精度。

1 二維ANCF 橫向高階剪切梁單元

Li 等引入曲率坐標(biāo)的同時，共用插值項系數(shù)，引入了橫向高階插值項-3x2y2，建立了一種新的二維ANCF 橫向高階剪切梁單元[12]，文中簡記為單元LPF。本文在此基礎(chǔ)上再一次通過共用x3系數(shù)，改為引入不同橫向高階插值項y2，構(gòu)建了新的二維高階梁單元，記為E1，其位置矢量為：

根據(jù)文獻(xiàn)[10]，低階單元橫縱向插值不一致使得其并不滿足廣義胡克定律。

根據(jù)應(yīng)力應(yīng)變關(guān)系σ=Dε，有

不考慮橫縱向高階耦合項時有

聯(lián)立式（2）和式（3），有

根據(jù)橫縱向應(yīng)變的定義，有

因此，新單元中通過共用系數(shù)引入的橫向高階插值項y2并能夠滿足廣義胡克定律，可以較好地緩解泊松閉鎖問題，提高收斂精度。

2 動力學(xué)方程

絕對節(jié)點(diǎn)坐標(biāo)法描述的二維梁單元動力學(xué)方程形式如下：

式中：Me——單元質(zhì)量矩陣；Ke（e）——單元剛度矩陣；Qe——單元廣義外力列陣。

2.1 單元質(zhì)量矩陣

式中：S——形函數(shù)矩陣；e——節(jié)點(diǎn)廣義坐標(biāo)向量。質(zhì)量矩陣由單元動能推導(dǎo)得到

因此，單元質(zhì)量矩陣為

2.2 單元剛度矩陣

一般連續(xù)介質(zhì)力學(xué)方法是ANCF 方法的基本力學(xué)原理。由格林-拉格朗日應(yīng)變張量可得：

其中任意一點(diǎn)變形梯度為

由2 階張量的voigt 標(biāo)記，將應(yīng)變張量改寫為應(yīng)變向量形式：

這里，泊松比v 會使得橫縱向正應(yīng)變ε11和ε22發(fā)生耦合。

彈性力列陣：

式中：K（e）——剛度矩陣。

2.3 廣義外力列陣

單元在受到外力F 作用下，單元廣義外力為

3 數(shù)值實例

本節(jié)通過懸臂梁靜力學(xué)、動力學(xué)等數(shù)值實例，對新單元的力學(xué)性能進(jìn)行測試。

3.1 靜力學(xué)大變形

懸臂梁如圖1 所示，尺寸為：長l=2 m，寬w=0.1 m，高h(yuǎn)=0.1 m。楊氏模量E=2.07×1011Pa，泊松比v=0.3，在懸臂端點(diǎn)重力方向施加的集中力F=500 000 N。以有限元軟件ANSYS 中Beam 188單元仿真結(jié)果為參照。

圖1 大變形懸臂梁Fig.1 Large deformation cantilever beam

表1 是不同單元在不同單元數(shù)劃分下，懸臂梁端點(diǎn)處y 方向位移情況，圖2 是不同單元與ANSYS 計算結(jié)果的誤差情況。由表1 和圖2 可知，新單元E1 與原單元LPF 相比，收斂精度有明顯提高，泊松閉鎖問題得到有效緩解。單元數(shù)增加時，新單元與ANSYS 計算結(jié)果誤差接近于0。

表1 大變形下末端點(diǎn)Y 方向位移計算結(jié)果Tab.1 Calculation results of the displacement in Y direction of end point under large deformation

圖2 各單元與ANSYS 在Y 方向位移誤差Fig.2 Displacement error of each element in Y direction compared with ANSYS

3.3 單擺

單擺如圖3 所示，長l=1 m，寬w=0.006 32 m，高h(yuǎn)=0.252 98 m，密度ρ=5 540 kg/m3，泊松比v=0.3，楊氏模量E=700 000 Pa，重力g 為外部載荷。計算0～1.1 s 單擺的運(yùn)動情況，在MATLAB 平臺進(jìn)行仿真計算。

圖3 單擺模型Fig.3 Simple pendulum model

圖4 是不同時刻各單元位姿，圖5、圖6 是各單元末端點(diǎn)不同時刻下的x、y 方向位移情況。根據(jù)圖4 可知，不同時刻下，新單元的位姿與ABAQUS 軟件仿真結(jié)果基本一致。由圖5、圖6 可知，新單元E1 末端點(diǎn)不同時刻下X、Y 方向的位移情況與ABAQUS 末端點(diǎn)的位移也基本吻合。而單元LPF 與ABAQUS有限元軟件的結(jié)果存在較大誤差，故新單元能夠更精確地描述動力學(xué)大變形運(yùn)動過程。

圖4 不同時刻各單元位姿Fig.4 Poses of each element at different moments,

圖5 不同時刻各單元末端點(diǎn)X 方向位移Fig.5 Displacement of end points of each element in X direction at different time

圖6 不同時刻各單元末端點(diǎn)Y 方向位移Fig.6 Displacement of end points of each element in Y direction at different time

4 結(jié)論

本文提出了一種共用系數(shù)的二維ANCF 高階剪切梁單元，通過理論推導(dǎo)結(jié)合靜力學(xué)、動力學(xué)數(shù)值實例，證明了（1）新單元相較于單元LPF 能夠有效緩解泊松閉鎖問題，提高單元靜力學(xué)大變形中的的收斂精度。（2）相較于單元LPF，新單元能夠更加精確地進(jìn)行動力學(xué)仿真和分析。