郭嘉輝，趙春花，獨亞平，周 川，張立強

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

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

Shabana提出的絕對節(jié)點坐標法(absolute nodal coordinate formulation，ANCF)是柔性多體系統(tǒng)動力學建模方法之一,單元節(jié)點廣義坐標定義在全局坐標系下，使得單元質(zhì)量矩陣為常數(shù)陣，且單元運動微分方程組裝成系統(tǒng)總運動微分方程時，無需進行坐標轉(zhuǎn)換，為其動力學模型求解以及復雜系統(tǒng)動力學建模提供了極大便利[1-3]。另外，該方法用位移梯度代替?zhèn)鹘y(tǒng)有限單元中的有限轉(zhuǎn)動，并丟棄了各種假設(shè)，因此能模擬出柔性部件的真實變形，為處理大變形、大范圍運動問題提供了可能。絕對節(jié)點坐標剪切梁單元位移模式通常采用多項式作為近似函數(shù)，那么多項式位移模式的選擇是決定梁單元力學特性的關(guān)鍵。絕對節(jié)點坐標剪切梁單元首先提出的都是橫向一次剪切梁單元，這些單元[4-9]的位移模式在縱向分量上都是高次插值，在橫向分量上都是線性插值，因此又稱為低階梁單元。沈振興等[10]和趙春花等[11]分別針對三維低階梁單元和二維低階梁單元，理論分析了線彈性材料的泊松閉鎖問題，指出其主要原因是多項式位移模式橫向一次和縱向高次不滿足本構(gòu)方程表達的縱橫向應變之間關(guān)于泊松比的線性關(guān)系，導致低階梁單元求解結(jié)果精度低。2010年以來，絕對節(jié)點坐標高階梁單元即橫向高次的剪切梁單元得到研究[12-16]。目前，已開發(fā)的高階梁單元在原橫向一次的基礎(chǔ)上均只引入了橫向二次多項式。因為橫向二次多項式的存在，其位移模式具備了二次完全多項式，滿足了本構(gòu)方程表達的縱橫向應變之間關(guān)于泊松比的線性關(guān)系，它們的精度被證明明顯好于低階梁單元。

目前，橡膠材料本構(gòu)關(guān)系主要采用3種不可壓縮模型：Neo-Hookean，Mooney-Rivlin和Yeoh。Maquedahe等[17]、Jung等[18]、Orzechowski等[19-20]和Xu等[21-22]采用絕對節(jié)點坐標三維低階梁單元和三維高階梁單元模擬橡膠梁的彎曲大變形，計算精度較高，但由于自由度較多，且非線性材料本身收斂迭代次數(shù)較多，導致計算時間較長，對計算機硬件有較高要求。目前，涉及二維絕對節(jié)點坐標剪切梁單元在超彈性材料非線性方面的力學特性研究較少。二維絕對節(jié)點坐標剪切梁單元的自由度數(shù)明顯少于三維高階梁單元。在可以轉(zhuǎn)化為二維問題的情況下，使用二維絕對節(jié)點坐標剪切梁單元對提高計算效率具有很大益處。因此課題組基于絕對節(jié)點坐標高階梁單元研究了大變形橡膠梁力學模型的計算精度和計算效率，探索了3種不可壓縮超彈性模型(Neo-Hookean，Mooney-Rivlin和Yeoh)模擬橡膠梁大變形力學特性的有效性，以及單元自由度數(shù)對計算效率的影響。

1 二維絕對節(jié)點坐標橫向高階梁單元

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

表1 二維橫向高階ANCF梁單元Table 1 Two-dimensional transverse high-order ANCF beam element

(1)

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

拉格朗日方程是建立單元運動微分方程的常用方法，其表達式如下:

(2)

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

在利用拉格朗日方程建立單元運動微分方程時需要建立單元的動能T和應變能U微分方程,且有：

(3)

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

(4)

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

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

(5)

式中：QT為單元彈性力列陣，由單元應變能U對單元節(jié)點廣義坐標e求偏導得到。

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

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

2.1 Neo-Hookean模型

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

(6)

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

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

(7)

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

(8)

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

(9)

2.2 Mooney-Rivlin模型

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

(10)

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

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

(11)

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

(12)

2.3 Yeoh模型

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

(13)

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

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

(14)

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

(15)

3 算例分析

課題組利用靜態(tài)和動態(tài)大變形實例研究了3種不可壓縮超彈性模型Neo-Hookean，Mooney-Rivlin和Yeoh，以及二維絕對節(jié)點坐標高階梁單元模擬橡膠梁大變形力學特性的有效性。

3.1 靜態(tài)大變形

靜態(tài)大變形實例選用的是受重力作用的橡膠懸臂梁。靜態(tài)懸臂梁模型如圖1所示，橡膠梁左端固定，另一端自由，整個梁在自身重力作用下產(chǎn)生大變形。重力常數(shù)g=9.81 m/s2，橡膠梁長l=160 mm，寬w=7 mm，高h=5 mm。橡膠梁材料參數(shù)如表2所示。其中Neo-Hookean模型與Yeoh模型參數(shù)與文獻[18]中使用參數(shù)相同，材料密度為2 150 kg/m3，Mooney-Rivlin模型參數(shù)與文獻[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 基于不可壓縮超彈性模型的二維高階梁單元計算懸臂梁位移收斂速度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所示為各個單元在不可壓縮超彈性模型下的懸臂梁自由端部縱橫向位移的收斂值，結(jié)果與Omar提出的二維絕對節(jié)點坐標低階單元[4]566進行對比，表中用LB2n-12表示該單元模型，亦與ABAQUS的結(jié)果對比。從表3可以看出：①二維絕對節(jié)點坐標低階單元在3種不可壓縮超彈性模型下的結(jié)果與高階單元和ABAQUS軟件的結(jié)果相比偏小很多；②對于Neo-Hookean模型，二維絕對節(jié)點坐標高階單元的結(jié)果與ABAQUS軟件的結(jié)果相差較大；③對于Mooney-Rivlin模型，二維絕對節(jié)點坐標高階單元的結(jié)果與ABAQUS軟件的結(jié)果比較接近；④對于Yeoh模型，二維絕對節(jié)點坐標高階單元的結(jié)果與ABAQUS軟件的結(jié)果非常接近?？梢奩eoh模型是最適合超彈性大變形分析的，二維絕對節(jié)點坐標高階梁單元從收斂精度的角度最適合超彈性大變形分析。

表3 3種不可壓縮超彈性模型下二維高階梁單元計算懸臂梁位移精度Table 3 Convergence precision of two dimensional high-order beam elements with three incompressible hyperelastic models

3.2 動態(tài)大變形

動態(tài)大變形實例選用的是受重力作用的柔性單擺，具體模型如圖3所示。單擺的物理參數(shù)：長l=350 mm，寬w=7 mm，高h=5 mm。材料參數(shù)與靜力學相同，具體參數(shù)如表2所示。單擺的一端通過銷與地面連接。單擺初始位置水平，初速度為零。重力常數(shù)g=9.81 m/s2。本算例中的柔性單擺在擺動過程中會產(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個單元分析Neo-Hookean，Mooney-Rivlin和Yeoh的3種不可壓縮超彈性模型的計算精度，結(jié)果如圖5所示。從圖5可以看出：①對于Neo-Hookean和Mooney-Rivlin模型，二維絕對節(jié)點坐標高階單元的結(jié)果在0.2，0.3，0.4和0.5 s，變形越來越大時與ABAQUS軟件的結(jié)果相差越來越大；②對于Yeoh模型，二維絕對節(jié)點坐標高階單元的結(jié)果在0.0，0.1，0.2，0.3，0.4和0.5 s，與ABAQUS軟件的結(jié)果都非常接近。因此，Yeoh模型是最適合超彈性大變形分析的。從圖5(c)可以看出，二維絕對節(jié)點坐標高階梁單元從收斂精度的角度都適合超彈性大變形分析。

圖5 基于3種不可壓縮超彈性模型各單元模型模擬的橡膠單擺在不同時刻的構(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 計算效率分析

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

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

表4 各單元CPU運行時間Table 4 CPU running time of each element

4 結(jié)論

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

1)Yeoh模型是最適合超彈性大變形分析；從計算精度角度二維絕對節(jié)點坐標高階梁單元都適合超彈性大變形分析，但它們的計算效率相差較大。

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

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

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