絕對節點坐標四面體單元建模與動力學分析

鮑康文， 趙春花， 陶羽玲

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

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

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

課題組對四面體單元進行了仔細研究，推導了動能和應變能的體積坐標表達式，并利用拉格朗日方程建立了運動微分方程。通過對靜態和動態實例進行分析，與ABAQUS軟件中的四面體單元進行對比，討論了絕對節點坐標四面體單元的收斂性和精度。

1 笛卡爾坐標和體積坐標

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

(1)

四面體單元有4個節點，4個節點的位置矢量表示為

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

式中：xK,yk和zk代表節點k定義在全局參考系下的笛卡爾坐標。

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

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

顯見，這4個體積坐標并非獨立，存在如下約束關系：

ξ+η+ζ+χ=1。

(2)

因此四面體單元任意點的笛卡爾位置坐標可用體積坐標表示為：

(3)

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

從式(3)可以看出，笛卡爾坐標和體積坐標之間是線性映射關系，笛卡爾位置坐標是體積坐標的函數。其中4個節點對應的體積坐標為(1,0,0,0)，(0,1,0,0)，(0,0,1,0)和(0,0,0,1)。

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

(4)

式中Δ=|V|/6。

其中：

(5)

(7)

(8)

2 絕對節點坐標四面體單元

絕對節點坐標4節點四面體單元每個節點有12個自由度，每個四面體單元有48個自由度。笛卡爾坐標系下的單元節點坐標列陣為

體積坐標下的單元節點坐標列陣用ek表示。但體積坐標下各個節點的坐標列陣不完全一樣，如圖2所示。

圖2 絕對節點坐標四面體單元模型Figure 2 Tetrahedral element model with ANCF

根據絕對節點坐標法的定義，四面體單元內任意一點的位置場可以表示為

r=Se。

(8)

式中：S是絕對節點坐標法4節點四面體單元的形函數矩陣，是體積坐標的函數，是3×48的矩陣;e是體積坐標系下的單元節點坐標矢量。

其中：

S=[s1Is2Is3Is4Is5Is6Is7Is8Is9Is10Is11Is12Is13Is14Is15Is16I];

(9)

(10)

式中I是3×3的單位矩陣，采用文獻[9]給出的形函數定義。

根據笛卡爾坐標系與體積坐標系下的位置矢量梯度關系，可得體積坐標系下的單元節點坐標矢量e和笛卡爾坐標系下的單元節點坐標矢量p的轉換矩陣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)

因為體積坐標下各個節點的坐標列陣不完全一樣，因此對應的轉換矩陣也不同，4個節點對應4個子轉換矩陣如下：

(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。

總轉換矩陣T是4個子轉換矩陣的組合，是48×48的矩陣：

(18)

3 體積坐標下的動力學模型

3.1 動能

根據式(8)和式(12)，可得體積坐標表示的動能：

(19)

因此單元質量矩陣

(20)

3.2 應變能

借助連續介質力學理論[12]構建絕對節點坐標法四面體單元的應變能。該單元的變形梯度矩陣用體積坐標表示，可以寫成如下形式：

(21)

其中：

(22)

式中：Siξ，Siη，Siζ和Siχ是形函數關于體積坐標的導數，i=1,2,3;Six,Siy和Siz是形函數關于笛卡爾坐標的導數，i=1,2,3。

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

(23)

其中：

(24)

根據線彈性材料的本構關系可得，應力矢量和應變矢量之間的關系為：

σ=Dε。

(25)

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

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

(26)

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

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

(27)

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

(28)

其中應變的體積坐標表達式可根據公式(23)得到。

3.3 單元運動微分方程

根據拉格朗日方程建立單元運動微分方程。將動能和應變能代入拉格朗日方程:

(29)

則單元運動微分方程為

(30)

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

(31)

其中：

廣義外力：

(32)

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

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

4 算例分析

4.1 懸臂梁靜力學分析

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

K(P)P=F。

(33)

由于該方程是節點坐標的非線性方程，筆者采用牛頓迭代法進行求解。

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

對于梁結構網格，采用Papplardo等[10]2921使用的網格策略，將5個四面體單元作為1個塊單元，采用圖4所示的2種塊單元網格。

圖4 四面體單元網格剖分圖Figure 4 Tetrahedral element mesh generation

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

圖5 大變形梁網格模型Figure 5 Beam mesh model in large deformation

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

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

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

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

4.2 轉動柔性梁動力學分析

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

圖7 轉動柔性梁模型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

仿真結果與ABAQUS軟件進行對比。如圖10所示，絕對節點坐標四面體單元顯示的精度更好，可以更有效地模擬柔性梁的大變形過程。

圖10 轉動柔性梁末端節點y向位移對比Figure 10 End displacement comparison of rotating flexible beam in y direction

5 結論

絕對節點坐標法對于大變形柔性多體系統建模在數值上已得到驗證。課題組基于拉格朗日方程，建立了絕對節點坐標四面體單元體積坐標表示的動力學模型。由于平面單元對復雜系統的適用性有局限，所以建立絕對節點坐標四面體單元模型為之后的柔性夾爪研究提供理論基礎。

課題組通過靜動態實例分析，將仿真結果與ABAQUS軟件中的一次四面體單元和二次四面體單元進行了對比分析。比較結果顯示，ANCF四面體單元的精度高于ABAQUS四面體一次單元，隨著單元數增加，精度可與二次四面體單元一致。證明了絕對節點坐標四面體單元的正確性。

隨著科學技術的發展，諸如橡膠等材料在柔性多體系統中得到廣泛應用，仍然采用線彈性本構模型的數值仿真不能滿足應用的需求，超彈性材料的應用有待進一步研究。