基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法的超彈性材料二維梁?jiǎn)卧獎(jiǎng)恿W(xué)研究

郭嘉輝，趙春花，獨(dú)亞平，周 川，張立強(qiáng)

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

隨著科學(xué)技術(shù)的發(fā)展，橡膠等超彈性材料正逐漸應(yīng)用于大變形柔性多體系統(tǒng)中，如輪胎、橡膠履帶、橡膠輸送帶、軟體夾持器和軟體機(jī)器人等。橡膠作為一種典型的超彈性材料，其材料特性表現(xiàn)為變形能力大、回彈快、無遲滯性以及體積近似不可壓縮等。超彈性材料具有非線性的應(yīng)力與應(yīng)變關(guān)系，相比傳統(tǒng)的線彈性材料，同時(shí)會(huì)產(chǎn)生幾何非線性和材料非線性的雙重復(fù)雜特性，給大變形柔性多體系統(tǒng)的動(dòng)力學(xué)仿真和力學(xué)特性分析提出了挑戰(zhàn)。

Shabana提出的絕對(duì)節(jié)點(diǎn)坐標(biāo)法(absolute nodal coordinate formulation，ANCF)是柔性多體系統(tǒng)動(dòng)力學(xué)建模方法之一,單元節(jié)點(diǎn)廣義坐標(biāo)定義在全局坐標(biāo)系下，使得單元質(zhì)量矩陣為常數(shù)陣，且單元運(yùn)動(dòng)微分方程組裝成系統(tǒng)總運(yùn)動(dòng)微分方程時(shí)，無需進(jìn)行坐標(biāo)轉(zhuǎn)換，為其動(dòng)力學(xué)模型求解以及復(fù)雜系統(tǒng)動(dòng)力學(xué)建模提供了極大便利[1-3]。另外，該方法用位移梯度代替?zhèn)鹘y(tǒng)有限單元中的有限轉(zhuǎn)動(dòng)，并丟棄了各種假設(shè)，因此能模擬出柔性部件的真實(shí)變形，為處理大變形、大范圍運(yùn)動(dòng)問題提供了可能。絕對(duì)節(jié)點(diǎn)坐標(biāo)剪切梁?jiǎn)卧灰颇Ｊ酵ǔ２捎枚囗?xiàng)式作為近似函數(shù)，那么多項(xiàng)式位移模式的選擇是決定梁?jiǎn)卧W(xué)特性的關(guān)鍵。絕對(duì)節(jié)點(diǎn)坐標(biāo)剪切梁?jiǎn)卧紫忍岢龅亩际菣M向一次剪切梁?jiǎn)卧?，這些單元[4-9]的位移模式在縱向分量上都是高次插值，在橫向分量上都是線性插值，因此又稱為低階梁?jiǎn)卧?。沈振興等[10]和趙春花等[11]分別針對(duì)三維低階梁?jiǎn)卧投S低階梁?jiǎn)卧?，理論分析了線彈性材料的泊松閉鎖問題，指出其主要原因是多項(xiàng)式位移模式橫向一次和縱向高次不滿足本構(gòu)方程表達(dá)的縱橫向應(yīng)變之間關(guān)于泊松比的線性關(guān)系，導(dǎo)致低階梁?jiǎn)卧蠼饨Y(jié)果精度低。2010年以來，絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧礄M向高次的剪切梁?jiǎn)卧玫窖芯縖12-16]。目前，已開發(fā)的高階梁?jiǎn)卧谠瓩M向一次的基礎(chǔ)上均只引入了橫向二次多項(xiàng)式。因?yàn)闄M向二次多項(xiàng)式的存在，其位移模式具備了二次完全多項(xiàng)式，滿足了本構(gòu)方程表達(dá)的縱橫向應(yīng)變之間關(guān)于泊松比的線性關(guān)系，它們的精度被證明明顯好于低階梁?jiǎn)卧?/p>

目前，橡膠材料本構(gòu)關(guān)系主要采用3種不可壓縮模型：Neo-Hookean，Mooney-Rivlin和Yeoh。Maquedahe等[17]、Jung等[18]、Orzechowski等[19-20]和Xu等[21-22]采用絕對(duì)節(jié)點(diǎn)坐標(biāo)三維低階梁?jiǎn)卧腿S高階梁?jiǎn)卧M橡膠梁的彎曲大變形，計(jì)算精度較高，但由于自由度較多，且非線性材料本身收斂迭代次數(shù)較多，導(dǎo)致計(jì)算時(shí)間較長，對(duì)計(jì)算機(jī)硬件有較高要求。目前，涉及二維絕對(duì)節(jié)點(diǎn)坐標(biāo)剪切梁?jiǎn)卧诔瑥椥圆牧戏蔷€性方面的力學(xué)特性研究較少。二維絕對(duì)節(jié)點(diǎn)坐標(biāo)剪切梁?jiǎn)卧淖杂啥葦?shù)明顯少于三維高階梁?jiǎn)卧?。在可以轉(zhuǎn)化為二維問題的情況下，使用二維絕對(duì)節(jié)點(diǎn)坐標(biāo)剪切梁?jiǎn)卧獙?duì)提高計(jì)算效率具有很大益處。因此課題組基于絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧芯苛舜笞冃蜗鹉z梁力學(xué)模型的計(jì)算精度和計(jì)算效率，探索了3種不可壓縮超彈性模型(Neo-Hookean，Mooney-Rivlin和Yeoh)模擬橡膠梁大變形力學(xué)特性的有效性，以及單元自由度數(shù)對(duì)計(jì)算效率的影響。

1 二維絕對(duì)節(jié)點(diǎn)坐標(biāo)橫向高階梁?jiǎn)卧?/h2>

絕對(duì)節(jié)點(diǎn)坐標(biāo)法在絕對(duì)坐標(biāo)系下假設(shè)單元內(nèi)任意一點(diǎn)位置坐標(biāo)的多項(xiàng)式位移場(chǎng)。目前存在3種二維絕對(duì)節(jié)點(diǎn)坐標(biāo)橫向高階梁?jiǎn)卧?，分別由Patel等[16]2929和Zhao等[11]482提出。它們的多項(xiàng)式位移模式如表1所示。其中：ri為單元內(nèi)任意一點(diǎn)的位置坐標(biāo)，二維情況下，i=1,2；ai1～ai8為多項(xiàng)式位移模式的未知系數(shù)。利用節(jié)點(diǎn)廣義坐標(biāo)，推導(dǎo)得到單元形函數(shù)，這樣任意一點(diǎn)的位置坐標(biāo)矢量r可表示成形函數(shù)矩陣S和單元廣義坐標(biāo)列陣e的表達(dá)式:

表1 二維橫向高階ANCF梁?jiǎn)卧猅able 1 Two-dimensional transverse high-order ANCF beam element

(1)

式中：形函數(shù)S是單元坐標(biāo)x和y的函數(shù)，是2×n矩陣函數(shù)，n代表單元自由度數(shù)。

拉格朗日方程是建立單元運(yùn)動(dòng)微分方程的常用方法，其表達(dá)式如下:

(2)

式中：QK為廣義外力列陣;t是時(shí)間。

在利用拉格朗日方程建立單元運(yùn)動(dòng)微分方程時(shí)需要建立單元的動(dòng)能T和應(yīng)變能U微分方程,且有：

(3)

根據(jù)公式(3)得單元的質(zhì)量矩陣為：

(4)

式中ST是形函數(shù)矩陣的轉(zhuǎn)置。

超彈性材料單元應(yīng)變能U是由單元應(yīng)變能密度函數(shù)積分得到。將動(dòng)能和應(yīng)變能代入拉格朗日方程，則單元的運(yùn)動(dòng)微分方程為：

(5)

式中：QT為單元彈性力列陣，由單元應(yīng)變能U對(duì)單元節(jié)點(diǎn)廣義坐標(biāo)e求偏導(dǎo)得到。

2 不可壓縮超彈性本構(gòu)彈性力模型

橡膠材料的應(yīng)力與應(yīng)變關(guān)系是非線性的，其材料力學(xué)特性目前主要由應(yīng)變能密度函數(shù)來表示[23]。通過對(duì)應(yīng)變能密度函數(shù)積分得到應(yīng)變能，應(yīng)變能對(duì)單元廣義坐標(biāo)列陣e求偏導(dǎo)得到彈性力列陣。為了保證單元的不可壓縮性，課題組采用了罰函數(shù)法。

2.1 Neo-Hookean模型

不可壓縮Neo-Hookean模型的應(yīng)變能密度函數(shù)為:

(6)

式中：C=JTJ是右格林柯西應(yīng)變張量;J是位置梯度矩陣;tr(C)是矩陣C的跡;det (J)是矩陣J的行列式，為了滿足單元的不可壓縮性，通過不斷的調(diào)整懲罰系數(shù)α，使得det (J)=1;C10是材料系數(shù)。

對(duì)于初始是直梁的情況，根據(jù)方程(1)，單元內(nèi)任意一點(diǎn)的位置梯度：

(7)

式中：r1x和r1y是r1分別對(duì)x和y的導(dǎo)數(shù)；r2x和r2y是r2分別對(duì)x和y的導(dǎo)數(shù)；S1x和S1y是S1分別對(duì)x和y的導(dǎo)數(shù)；S2x和S2y是S2分別對(duì)x和y的導(dǎo)數(shù)。將式(7)代入式(6)，并進(jìn)行積分得不可壓縮Neo-Hookean模型的應(yīng)變能為：

(8)

將方程(8)對(duì)單元節(jié)點(diǎn)廣義坐標(biāo)列陣e求導(dǎo)，得單元彈性力列陣QT，并根據(jù)方程(7)可進(jìn)一步將單元彈性力列陣表示為：

(9)

2.2 Mooney-Rivlin模型

不可壓縮Mooney-Rivlin模型的應(yīng)變能密度函數(shù)為：

(10)

式中C10，C01是材料系數(shù)。

將方程(7)代入方程(10)并進(jìn)行積分，得不可壓縮Mooney-Rivlin模型的應(yīng)變能為：

(11)

將方程(11)對(duì)單元節(jié)點(diǎn)廣義坐標(biāo)列陣e求導(dǎo)，并根據(jù)方程(7)進(jìn)一步推導(dǎo)出單元彈性力列陣：

(12)

2.3 Yeoh模型

不可壓縮Yeoh模型的應(yīng)變能密度函數(shù)為：

(13)

式中C10，C20和C30是材料系數(shù)。

將方程(7)代入方程(13)，并進(jìn)行積分得不可壓縮Yeoh模型的應(yīng)變能為：

(14)

將方程(14)對(duì)單元節(jié)點(diǎn)廣義坐標(biāo)列陣e求導(dǎo)，并根據(jù)方程(7)進(jìn)一步推導(dǎo)單元彈性力列陣：

(15)

3 算例分析

課題組利用靜態(tài)和動(dòng)態(tài)大變形實(shí)例研究了3種不可壓縮超彈性模型Neo-Hookean，Mooney-Rivlin和Yeoh，以及二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧M橡膠梁大變形力學(xué)特性的有效性。

3.1 靜態(tài)大變形

靜態(tài)大變形實(shí)例選用的是受重力作用的橡膠懸臂梁。靜態(tài)懸臂梁模型如圖1所示，橡膠梁左端固定，另一端自由，整個(gè)梁在自身重力作用下產(chǎn)生大變形。重力常數(shù)g=9.81 m/s2，橡膠梁長l=160 mm，寬w=7 mm，高h(yuǎn)=5 mm。橡膠梁材料參數(shù)如表2所示。其中Neo-Hookean模型與Yeoh模型參數(shù)與文獻(xiàn)[18]中使用參數(shù)相同，材料密度為2 150 kg/m3，Mooney-Rivlin模型參數(shù)與文獻(xiàn)[19]中使用參數(shù)相同，材料密度為7 200 kg/m3，選取懲罰系數(shù)α=1 000 MPa。

圖1 靜態(tài)懸臂梁模型Figure 1 Static cantilever beam model

表2 各模型的材料參數(shù)Table 2 Material parameters of each model

3.1.1 收斂速度

在采用Neo-Hookean，Mooney-Rivlin和Yeoh的3種不可壓縮超彈性模型下，HB2n-12，HB3n-12和HB2n-16的收斂情況如圖2所示。從圖2可以看出，在Neo-Hookean，Mooney-Rivlin和Yeoh的3種不可壓縮超彈性模型下，HB3n-12收斂需要的單元數(shù)最少，HB2n-12次之，HB2n-16最多。

圖2 基于不可壓縮超彈性模型的二維高階梁?jiǎn)卧?jì)算懸臂梁位移收斂速度Figure 2 Convergence rates of cantilever beam displacement calculated by two-dimensional high-order beam element based on incompressible hyperelastic model

3.1.2 收斂精度

表3所示為各個(gè)單元在不可壓縮超彈性模型下的懸臂梁自由端部縱橫向位移的收斂值，結(jié)果與Omar提出的二維絕對(duì)節(jié)點(diǎn)坐標(biāo)低階單元[4]566進(jìn)行對(duì)比，表中用LB2n-12表示該單元模型，亦與ABAQUS的結(jié)果對(duì)比。從表3可以看出：①二維絕對(duì)節(jié)點(diǎn)坐標(biāo)低階單元在3種不可壓縮超彈性模型下的結(jié)果與高階單元和ABAQUS軟件的結(jié)果相比偏小很多；②對(duì)于Neo-Hookean模型，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階單元的結(jié)果與ABAQUS軟件的結(jié)果相差較大；③對(duì)于Mooney-Rivlin模型，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階單元的結(jié)果與ABAQUS軟件的結(jié)果比較接近；④對(duì)于Yeoh模型，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階單元的結(jié)果與ABAQUS軟件的結(jié)果非常接近?？梢奩eoh模型是最適合超彈性大變形分析的，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧獜氖諗烤鹊慕嵌茸钸m合超彈性大變形分析。

表3 3種不可壓縮超彈性模型下二維高階梁?jiǎn)卧?jì)算懸臂梁位移精度Table 3 Convergence precision of two dimensional high-order beam elements with three incompressible hyperelastic models

3.2 動(dòng)態(tài)大變形

動(dòng)態(tài)大變形實(shí)例選用的是受重力作用的柔性單擺，具體模型如圖3所示。單擺的物理參數(shù)：長l=350 mm，寬w=7 mm，高h(yuǎn)=5 mm。材料參數(shù)與靜力學(xué)相同，具體參數(shù)如表2所示。單擺的一端通過銷與地面連接。單擺初始位置水平，初速度為零。重力常數(shù)g=9.81 m/s2。本算例中的柔性單擺在擺動(dòng)過程中會(huì)產(chǎn)生較大變形。

圖3 柔性單擺模型Figure 3 Flexible pendulum model

3.2.1 收斂速度

在采用Neo-Hookean，Mooney-Rivlin和Yeoh的3種不可壓縮超彈性模型下，HB2n-12，HB3n-12和HB2n-16的收斂情況是基本一致，所以只展示了圖4所示的Yeoh不可壓縮超彈性模型下，HB2n-12，HB3n-12和HB2n-16的收斂情況。從圖4可以看出，HB3n-12收斂需要的單元數(shù)是16，而HB2n-12和 HB2n-16收斂需要的單元數(shù)是32。

圖4 基于Yeoh本構(gòu)方程不同單元模型的收斂速度Figure 4 Convergence rates of different element models based on Yeoh constitutive equation

3.2.2 收斂精度

根據(jù)HB2n-12，HB3n-12和HB2n-16收斂需要的單元數(shù)，選擇了32個(gè)單元分析Neo-Hookean，Mooney-Rivlin和Yeoh的3種不可壓縮超彈性模型的計(jì)算精度，結(jié)果如圖5所示。從圖5可以看出：①對(duì)于Neo-Hookean和Mooney-Rivlin模型，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階單元的結(jié)果在0.2，0.3，0.4和0.5 s，變形越來越大時(shí)與ABAQUS軟件的結(jié)果相差越來越大；②對(duì)于Yeoh模型，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階單元的結(jié)果在0.0，0.1，0.2，0.3，0.4和0.5 s，與ABAQUS軟件的結(jié)果都非常接近。因此，Yeoh模型是最適合超彈性大變形分析的。從圖5(c)可以看出，二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧獜氖諗烤鹊慕嵌榷歼m合超彈性大變形分析。

圖5 基于3種不可壓縮超彈性模型各單元模型模擬的橡膠單擺在不同時(shí)刻的構(gòu)型圖Figure 5 Configuration diagrams of rubber pendulum at different moments simulated by each element model based on three incompressible hyperelastic models

3.3 計(jì)算效率分析

從3.1和3.2節(jié)可知，Yeoh模型是最適合橡膠梁大變形力學(xué)特性分析的。橡膠梁大變形力學(xué)特性分析需要考慮梁的幾何非線性和材料非線性，因此對(duì)模型的計(jì)算效率提出了比較大的挑戰(zhàn)。已知3種二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧諗克鑶卧獢?shù)是不同的，HB3n-12收斂所需要單元數(shù)最少，最容易收斂，HB2n-12次之，HB2n-16收斂所需要單元數(shù)最多。這與單元的構(gòu)造即單元位移模式和節(jié)點(diǎn)廣義坐標(biāo)類型有關(guān)。

此外，相同單元數(shù)下，3個(gè)單元的CPU運(yùn)行時(shí)間也是不一樣的。表4所示為HB2n-12，HB3n-12和HB2n-16計(jì)算動(dòng)態(tài)大變形實(shí)例0.5 s所消耗的CPU運(yùn)行時(shí)間(仿真環(huán)境為8 GB運(yùn)行內(nèi)存，Inter Core i5-5200U，CPU主頻為2.20 GHz的計(jì)算平臺(tái))。從表4可以看出，在單元數(shù)分別是4，8，16和32時(shí)，HB3n-12和HB2n-12所消耗CPU運(yùn)行時(shí)間相差不大，HB3n-12比HB2n-12多一些時(shí)間，而HB2n-16所消耗CPU運(yùn)行時(shí)間就遠(yuǎn)多于HB3n-12和HB2n-12，雖然HB2n-16的自由度數(shù)16只比HB3n-12和HB2n-12的自由度數(shù)12多了4，但CPU運(yùn)行時(shí)間已經(jīng)達(dá)到它們的2倍左右。因此單元自由度數(shù)對(duì)計(jì)算效率影響較大。

表4 各單元CPU運(yùn)行時(shí)間Table 4 CPU running time of each element

4 結(jié)論

課題組基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法研究了HB2n-12，HB3n-12和HB2n-16的3種二維高階單元和Neo-Hookean，Mooney-Rivlin和Yeoh的3種不可壓縮超彈性模型模擬橡膠梁大變形的力學(xué)特性?；陟o態(tài)和動(dòng)態(tài)大變形實(shí)例，分析了3種單元的計(jì)算精度和計(jì)算效率以及3種不可壓縮超彈性模型模擬橡膠梁大變形力學(xué)特性的有效性。具體結(jié)論如下：

1)Yeoh模型是最適合超彈性大變形分析；從計(jì)算精度角度二維絕對(duì)節(jié)點(diǎn)坐標(biāo)高階梁?jiǎn)卧歼m合超彈性大變形分析，但它們的計(jì)算效率相差較大。

2)HB3n-12收斂需要的單元數(shù)最少，最容易收斂，HB2n-12次之，HB2n-16收斂所需要單元數(shù)最多。這與單元的構(gòu)造即單元位移模式和節(jié)點(diǎn)廣義坐標(biāo)類型有關(guān)。

3)相同單元數(shù)下，HB2n-16所消耗CPU運(yùn)行時(shí)間是HB3n-12和HB2n-12所消耗的2倍左右，因此單元自由度數(shù)對(duì)CPU運(yùn)行時(shí)間影響較大。

因此，橡膠梁大變形分析適合采用Yeoh模型，單元多項(xiàng)式位移模式、節(jié)點(diǎn)廣義坐標(biāo)類型以及單元自由度數(shù)對(duì)單元計(jì)算效率影響較大。