鮑康文， 趙春花， 陶羽玲

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

近年來，隨著科學(xué)技術(shù)的快速發(fā)展，柔性構(gòu)件得到了廣泛應(yīng)用，柔性多體系統(tǒng)的研究也受到了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注[1]。對(duì)于具有復(fù)雜外形、不能簡(jiǎn)化為梁板殼的不規(guī)則柔性構(gòu)件，基于有限元理論的數(shù)值模擬只能采用體單元，而四面體單元[2]是一種常用的體單元。絕對(duì)節(jié)點(diǎn)坐標(biāo)法(ANCF)由Shabana等于1996年提出[3-4]，是用于描述柔性多體系統(tǒng)動(dòng)力學(xué)特性的一種建模方法。近20年來，絕對(duì)節(jié)點(diǎn)坐標(biāo)法已在數(shù)值和實(shí)驗(yàn)上得到驗(yàn)證，并成功地用于多種柔性多體系統(tǒng)的建模[5-6]。

但由于四面體單元的建模過程比較復(fù)雜，單元節(jié)點(diǎn)自由度較多，因此絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元的建模與應(yīng)用有待進(jìn)一步研究。Lan和Olshevskiy等[7-8]提出了一種ANCF實(shí)體四面體有限 元。該單元適合于復(fù)雜形狀結(jié)構(gòu)的建模，且自由度小，因而具有良好的性能，在結(jié)構(gòu)動(dòng)力學(xué)問題中有著廣泛的應(yīng)用。Mohamed[9]提出了一種完全參數(shù)化的絕對(duì)節(jié)點(diǎn)坐標(biāo)法四面體單元，該論文中的四面體單元形函數(shù)是采用笛卡爾坐標(biāo)構(gòu)造的。Pappalardo等[10]提出了ANCF等參四面體有限元的開發(fā)方法，討論了傳統(tǒng)有限元四面體單元坐標(biāo)參數(shù)化與ANCF四面體單元坐標(biāo)參數(shù)化的基本區(qū)別，并定義了表示體積參數(shù)和笛卡爾參數(shù)之間關(guān)系的單元轉(zhuǎn)換矩陣，從而可以使用標(biāo)準(zhǔn)的有限元裝配過程。Wang等[11]提出了一種不完全三次絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元，論文中還提出了采用面中心點(diǎn)和材料點(diǎn)的位置約束來構(gòu)建不完全三次單元。

課題組對(duì)四面體單元進(jìn)行了仔細(xì)研究，推導(dǎo)了動(dòng)能和應(yīng)變能的體積坐標(biāo)表達(dá)式，并利用拉格朗日方程建立了運(yùn)動(dòng)微分方程。通過對(duì)靜態(tài)和動(dòng)態(tài)實(shí)例進(jìn)行分析，與ABAQUS軟件中的四面體單元進(jìn)行對(duì)比，討論了絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元的收斂性和精度。

1 笛卡爾坐標(biāo)和體積坐標(biāo)

如圖1所示，四面體內(nèi)任意一點(diǎn)的位置場(chǎng)r=[r1r2r3]可以表示成笛卡爾參數(shù)x=[xyz]T的函數(shù)，如式(1)所示：

(1)

四面體單元有4個(gè)節(jié)點(diǎn)，4個(gè)節(jié)點(diǎn)的位置矢量表示為

vk=[xkykzk]T，k=1，2，3，4。

式中：xK,yk和zk代表節(jié)點(diǎn)k定義在全局參考系下的笛卡爾坐標(biāo)。

為了便于積分，四面體單元的笛卡爾坐標(biāo)需要轉(zhuǎn)換成體積坐標(biāo)[10]2910。四面體單元的4個(gè)節(jié)點(diǎn)分別用A,B,C和D表示，按照右螺旋規(guī)則排列，如圖1所示。四面體內(nèi)的任意點(diǎn)P將四面體劃分為4個(gè)子四面體PBCD，PCDA，PDAB和PABC。Δ代表四面體單元的體積,Δ1,Δ2,Δ3和Δ4代表4個(gè)子四面體的體積。則四面體單元的體積坐標(biāo)定義為

ξ=Δ1/Δ,η=Δ2/Δ,ζ=Δ3/Δ,χ=Δ4/Δ。

顯見，這4個(gè)體積坐標(biāo)并非獨(dú)立，存在如下約束關(guān)系：

ξ+η+ζ+χ=1。

(2)

因此四面體單元任意點(diǎn)的笛卡爾位置坐標(biāo)可用體積坐標(biāo)表示為：

(3)

圖1 四面體單元幾何模型Figure 1 Geometric model of tetrahedral element

從式(3)可以看出，笛卡爾坐標(biāo)和體積坐標(biāo)之間是線性映射關(guān)系，笛卡爾位置坐標(biāo)是體積坐標(biāo)的函數(shù)。其中4個(gè)節(jié)點(diǎn)對(duì)應(yīng)的體積坐標(biāo)為(1,0,0,0)，(0,1,0,0)，(0,0,1,0)和(0,0,0,1)。

根據(jù)式(2)和式(3)，體積坐標(biāo)也可表示為笛卡爾坐標(biāo)的函數(shù)：

(4)

式中Δ=|V|/6。

其中：

(5)

(7)

(8)

2 絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元

絕對(duì)節(jié)點(diǎn)坐標(biāo)4節(jié)點(diǎn)四面體單元每個(gè)節(jié)點(diǎn)有12個(gè)自由度，每個(gè)四面體單元有48個(gè)自由度。笛卡爾坐標(biāo)系下的單元節(jié)點(diǎn)坐標(biāo)列陣為

體積坐標(biāo)下的單元節(jié)點(diǎn)坐標(biāo)列陣用ek表示。但體積坐標(biāo)下各個(gè)節(jié)點(diǎn)的坐標(biāo)列陣不完全一樣，如圖2所示。

圖2 絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元模型Figure 2 Tetrahedral element model with ANCF

根據(jù)絕對(duì)節(jié)點(diǎn)坐標(biāo)法的定義，四面體單元內(nèi)任意一點(diǎn)的位置場(chǎng)可以表示為

r=Se。

(8)

式中：S是絕對(duì)節(jié)點(diǎn)坐標(biāo)法4節(jié)點(diǎn)四面體單元的形函數(shù)矩陣，是體積坐標(biāo)的函數(shù)，是3×48的矩陣;e是體積坐標(biāo)系下的單元節(jié)點(diǎn)坐標(biāo)矢量。

其中：

S=[s1Is2Is3Is4Is5Is6Is7Is8Is9Is10Is11Is12Is13Is14Is15Is16I];

(9)

(10)

式中I是3×3的單位矩陣，采用文獻(xiàn)[9]給出的形函數(shù)定義。

根據(jù)笛卡爾坐標(biāo)系與體積坐標(biāo)系下的位置矢量梯度關(guān)系，可得體積坐標(biāo)系下的單元節(jié)點(diǎn)坐標(biāo)矢量e和笛卡爾坐標(biāo)系下的單元節(jié)點(diǎn)坐標(biāo)矢量p的轉(zhuǎn)換矩陣T，形式如下:

e=Tp。

(11)

其中e可以表示為

e=[(e1)T(e2)T(e3)T(e4)T]。

(12)

p可以表示為

p=[(p1)T(p2)T(p3)T(p4)T]。

(13)

因?yàn)轶w積坐標(biāo)下各個(gè)節(jié)點(diǎn)的坐標(biāo)列陣不完全一樣，因此對(duì)應(yīng)的轉(zhuǎn)換矩陣也不同，4個(gè)節(jié)點(diǎn)對(duì)應(yīng)4個(gè)子轉(zhuǎn)換矩陣如下：

(14)

(15)

(16)

(17)

其中：

ai,j=xi-xj;bi,j=yi-yj;ci,j=zi-zj;

i=1,2,3,4,j=1,2,3,4,i≠j。

總轉(zhuǎn)換矩陣T是4個(gè)子轉(zhuǎn)換矩陣的組合，是48×48的矩陣：

(18)

3 體積坐標(biāo)下的動(dòng)力學(xué)模型

3.1 動(dòng)能

根據(jù)式(8)和式(12)，可得體積坐標(biāo)表示的動(dòng)能：

(19)

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

(20)

3.2 應(yīng)變能

借助連續(xù)介質(zhì)力學(xué)理論[12]構(gòu)建絕對(duì)節(jié)點(diǎn)坐標(biāo)法四面體單元的應(yīng)變能。該單元的變形梯度矩陣用體積坐標(biāo)表示，可以寫成如下形式：

(21)

其中：

(22)

式中：Siξ，Siη，Siζ和Siχ是形函數(shù)關(guān)于體積坐標(biāo)的導(dǎo)數(shù)，i=1,2,3;Six,Siy和Siz是形函數(shù)關(guān)于笛卡爾坐標(biāo)的導(dǎo)數(shù)，i=1,2,3。

根據(jù)式(21)四面體單元的變形梯度J，四面體單元的應(yīng)變張量采用拉格朗日格林應(yīng)變張量，可寫成如下形式：

(23)

其中：

(24)

根據(jù)線彈性材料的本構(gòu)關(guān)系可得，應(yīng)力矢量和應(yīng)變矢量之間的關(guān)系為：

σ=Dε。

(25)

式中D是材料的彈性模量矩陣。

對(duì)各向同性材料，矩陣D可以使用拉梅常數(shù)λ和μ來表示:

(26)

式中：ν為泊松比；E為彈性模量。

則應(yīng)變能可由彈性模量矩陣和拉格朗日應(yīng)變張量表示：

(27)

將式(23)和式(26)代入式(27)，可得體積坐標(biāo)表示的應(yīng)變能：

(28)

其中應(yīng)變的體積坐標(biāo)表達(dá)式可根據(jù)公式(23)得到。

3.3 單元運(yùn)動(dòng)微分方程

根據(jù)拉格朗日方程建立單元運(yùn)動(dòng)微分方程。將動(dòng)能和應(yīng)變能代入拉格朗日方程:

(29)

則單元運(yùn)動(dòng)微分方程為

(30)

Qk=Ke(p)p為單元彈性力矩陣，表示為：

(31)

其中：

廣義外力：

(32)

其中fe=[0 -ρVg0]T為單元重力。

式中：ρ為材料密度；V為四面體單元體積；g為重力加速度。

4 算例分析

4.1 懸臂梁靜力學(xué)分析

課題組采用懸臂梁結(jié)構(gòu)模型的長(zhǎng)度L=2 m，寬度W=0.1 m，高度H=0.1 m，末端節(jié)點(diǎn)施加的力F=500 kN，如圖3所示。模型的彈性模量E=207 GPa,泊松比ν=0.3。由于本算例屬于靜態(tài)非線性問題，不需要考慮單元的質(zhì)量矩陣，從而總運(yùn)動(dòng)微分方程表示為

K(P)P=F。

(33)

由于該方程是節(jié)點(diǎn)坐標(biāo)的非線性方程，筆者采用牛頓迭代法進(jìn)行求解。

圖3 懸臂梁模型Figure 3 Cantilever beam model

對(duì)于梁結(jié)構(gòu)網(wǎng)格，采用Papplardo等[10]2921使用的網(wǎng)格策略，將5個(gè)四面體單元作為1個(gè)塊單元，采用圖4所示的2種塊單元網(wǎng)格。

圖4 四面體單元網(wǎng)格剖分圖Figure 4 Tetrahedral element mesh generation

采用圖5所示大變形梁網(wǎng)格圖，該模型是由圖4所示大小相同的正方體塊組成。

圖5 大變形梁網(wǎng)格模型Figure 5 Beam mesh model in large deformation

表1為不同種類四面體單元收斂性分析。

表1 不同種類四面體單元收斂性分析

仿真結(jié)果與ABAQUS四面體一次單元、二次單元進(jìn)行對(duì)比，如表1和圖6所示。結(jié)果顯示，隨著單元數(shù)的增加，ANCF四面體單元的收斂性較好，并且，計(jì)算精度在單元數(shù)較少時(shí)候與ABAQUS一次四面體單元接近，但是當(dāng)單元數(shù)增加后，精度明顯提高，收斂結(jié)果與ABAQUS二次四面體單元的趨于一致。

圖6 懸臂梁末端位移Figure 6 End displacement of cantilever beam

4.2 轉(zhuǎn)動(dòng)柔性梁動(dòng)力學(xué)分析

在動(dòng)力學(xué)實(shí)例中，采用轉(zhuǎn)動(dòng)柔性梁，長(zhǎng)L=2 m，寬W=0.1 m，高H=0.1 m，如圖7所示。彈性模量E=700 kPa，泊松比ν=0.3，質(zhì)量密度ρ=5 540 kg/m3，重力加速度g=9.81 m/s2，仿真的時(shí)間為1.5 s，半個(gè)周期。如圖8～9所示，對(duì)絕對(duì)節(jié)點(diǎn)坐標(biāo)中單元數(shù)分別為200，400和600個(gè)四面體單元進(jìn)行動(dòng)力學(xué)收斂性分析。當(dāng)單元數(shù)為200時(shí)，模型在局部位置存在誤差。當(dāng)單元數(shù)增加到400和600時(shí)，模型已經(jīng)趨于收斂。

圖7 轉(zhuǎn)動(dòng)柔性梁模型Figure 7 Model of rotating flexible beam

圖8 ANCF四面體單元x方向收斂性分析Figure 8 Convergence analysis of ANCF tetrahedral element in x direction

圖9 ANCF四面體單元y方向收斂性分析Figure 9 Convergence analysis of ANCF tetrahedral element in y direction

仿真結(jié)果與ABAQUS軟件進(jìn)行對(duì)比。如圖10所示，絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元顯示的精度更好，可以更有效地模擬柔性梁的大變形過程。

圖10 轉(zhuǎn)動(dòng)柔性梁末端節(jié)點(diǎn)y向位移對(duì)比Figure 10 End displacement comparison of rotating flexible beam in y direction

5 結(jié)論

絕對(duì)節(jié)點(diǎn)坐標(biāo)法對(duì)于大變形柔性多體系統(tǒng)建模在數(shù)值上已得到驗(yàn)證。課題組基于拉格朗日方程，建立了絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元體積坐標(biāo)表示的動(dòng)力學(xué)模型。由于平面單元對(duì)復(fù)雜系統(tǒng)的適用性有局限，所以建立絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元模型為之后的柔性夾爪研究提供理論基礎(chǔ)。

課題組通過靜動(dòng)態(tài)實(shí)例分析，將仿真結(jié)果與ABAQUS軟件中的一次四面體單元和二次四面體單元進(jìn)行了對(duì)比分析。比較結(jié)果顯示，ANCF四面體單元的精度高于ABAQUS四面體一次單元，隨著單元數(shù)增加，精度可與二次四面體單元一致。證明了絕對(duì)節(jié)點(diǎn)坐標(biāo)四面體單元的正確性。

隨著科學(xué)技術(shù)的發(fā)展，諸如橡膠等材料在柔性多體系統(tǒng)中得到廣泛應(yīng)用，仍然采用線彈性本構(gòu)模型的數(shù)值仿真不能滿足應(yīng)用的需求，超彈性材料的應(yīng)用有待進(jìn)一步研究。