非協調六面體有限元分析方法研究

韋建平，韋宏法，韋志安，黃振之

(1.柳州孔輝汽車科技有限公司，廣西柳州 545006；2.柳州五菱柳機動力有限公司，廣西柳州 545005)

?

非協調六面體有限元分析方法研究

韋建平1，韋宏法1，韋志安2，黃振之1

(1.柳州孔輝汽車科技有限公司，廣西柳州 545006；2.柳州五菱柳機動力有限公司，廣西柳州 545005)

基于平面等參協調單元給出空間六面體有限元模型表達式。針對完備的協調單元由于剛度系數值比精確值大、導致在給定的載荷之下計算模型的變形比實際結構小的問題，建立了一套可以滿足分片檢驗條件的六面體非協調單元的有限元分析方法，有效解決了三維協調元為提高計算精度而導致的計算效率低、計算內存占用大的問題。并通過實例驗證了非協調元的計算精度。

等參元；非協調元；有限元；六面體網格

0 引言

協調等參元雖有良好的適應性和表達格式的簡明性，也因此得到廣泛的應用，但是從嚴格的意義上說，它的精度和效率仍是不夠高的。在三維坐標中，一次和二次完全多項式分別是4項和10項。而三維線性單元和二次單元卻分別具有8個和20個節點，也即三維等參元中有二分之一的節點自由度對計算精度是無貢獻的。因此，E L WILSON提出了非協調等參單元，對改進等參單元計算精度和提高計算效率是很有意義的[1-2]。

作者結合二維非協調有限元分析方法，并在文獻[3]提出的協調等參單元的基礎上，建立了一套可以滿足分片檢驗條件的三維六面體非協調單元的有限元解法，有效解決了三維非協調元的計算問題。

1 有限元力學基礎

有限元法中，需要用到彈性力學的基本方程和與之等效的變分原理。彈性體內任一點的位移可由沿直角坐標軸方向的3個位移分量u、v和w來表示，它的矩陣形式是[2]：

(1)

上式稱作位移列陣或位移向量。對于三維問題，彈性力學基本方程包括平衡方程、幾何方程(即應變和位移變化關系)、物理方程(即應力和應變關系)，還有力和幾何的邊界條件，在V域內，其方程一般形式如下：

平衡方程： Lσ+pv=0 (在V內)

幾何方程： ε=LTu (在V內)

物理方程： σ=Dε (在V內)

(2)

nσ=ps(在邊界V上)

2 協調六面體等參有限元

2.1 六面體單元概述

圖1(a)是一個物理空間中的任意四邊形單元。在該單元上建立一個參考坐標系ξη。四邊形單元的邊被橫坐標ξ和縱坐標η平分，它們所表示的方程是ξ=±1和η=±1，如圖1(a)所示。在該參考坐標系ξη下，該四邊形變換成一個邊長等于2的正方形，如圖1(b)所示。

圖1 四邊形等參數單元

同理，空間等參數單元可由平面問題等參單元推廣得到。空間八節點單元是直棱的，可以描述形狀復雜的三維結構體。八節點等參數單元及其參考坐標系ξηζ如圖2所示。在參考坐標系下，直棱的六面體是一個邊長等于2的立方體，坐標系ξηζ的原點位于形心處[1]。

圖2 八節點六面體等參單元

2.2 六面體單元特性

根據等參數單元的含義，它的位移和坐標都采用相同的形函數表示，即

位移函數：

(3)

坐標函數：

式中：ui、vi、wi(i=1,2,3,4,5,6,7,8)是整體坐標系xyz下節點的位移分量，xi、yi、zi是節點的整體坐標， Ni(ξ,η,ζ)是由參考坐標ξηζ表示的形函數。根據矩形單元的形函數公式，可以得出八節點等參單元的形函數Ni的表達式：

(4)

式中：ξi，ηi和ζi是節點i的局部坐標，分別位于六面體的8個角處，其對應坐標值分別為1或者-1，即：

(5)

圖2所示的單元節點中，節點編號1-2-3-4-5-6-7-8是逆時針并且由外向內順序的，而其參考坐標系ξηζ在整體坐標系xyz下的方向是由單元的節點編號決定的。

根據式(2)的幾何方程，等參六面體的應變ε公式為：

(6)

式(6)中，由于位移函數uvw和坐標函數xyz不在同一坐標系下，根據復合函數的求導法則，必須經過必要的轉換，記Ni,x、Ni,y、Ni,z分別表示形函數Ni對x、y和z的偏導數，它們與參考坐標系下對應的Ni,ξ、Ni,η、Ni,ζ關系式表示如下：

(7)

式中：矩陣J稱為雅克比(Jacobin)矩陣。其表達式如下所示：

根據式(6)、式(7)得出等參六面體的應變ε計算公式為：

(8)

式中：B稱為單元應變系數，從而單元應力的計算公式為：

(9)

式中：

(10)

其中：常數A1、A2和A3由式(11)定義：

(11)

綜上所述，并根據單元節點的剛度公式，三維等參單元的剛度矩陣表達式如下：

(12)

單元剛度矩陣Ke可劃分成n×n個子矩陣形式，其計算公式如下：

(i,j=1,2,…,8)

(13)

由于形函數對整體坐標x和y的導數中包含雅克比矩陣的逆矩陣，式(13)在一般情況下很難得到顯示方程，必須采用數值積分來進行求解，該方法可參考文獻[1]。

3 非協調元的凝聚

3.1 Wilson非協調元

對于三維線性單元來講，由于單元的插值函數中包含有非完全的ξη、ξζ和ηζ項，在用它表示純彎曲應力容易出現誤差。

而附加項的位移在三維線性單元的8個節點上都取零值時，它對節點位移沒有任何影響，而只對單元內部的位移起了調整作用。待定系數α1、α2和α3與單元邊界上的節點無關。

包含附加的無節點位移項的單元位移插值表示如下：

(14)

式中：Ni為式(4)所示的形函數。

等式(14)右邊由兩部分組成，前者為第2節介紹的具備協調性的等參單元位移形函數，后者為非協調部分的形函數，兩者矩陣維數不同，計算過程和計算方法是一致的。

因此，為了改善三維線性單元的性質，解決完備協調元中剛度系數值過大，提高其解精度，需將單元的非協調位移插值函數附加到內部無節點的位移項。

3.2 非協調元剛度的凝聚

根據上文所述，結合幾何關系和位移函數，加入非協調部分后，單元的應變矩陣公式變為如下式子：

(15)

將式(14)和式(15)等代入位移能泛函并按照通常步驟，此時可以得到位移、剛度與受力的矩陣式[3]：

(16)

其中：

(17)

式中：Pe為外載荷列陣。

對上述關系式進行展開求解，則單元內部位移、剛度矩陣以及載荷列陣之間的關系式

Keδe=Pe

(18)

式中：

此式即為包含附加內部位移項的單元剛度矩陣和載荷列陣。它是在原單元剛度矩陣和載荷列陣內增加了修正項而得到的。經過凝聚后，單元的自由度仍是原六面體單元的自由度，之后的分析計算步驟也跟之前的標準解題方法和步驟一樣。

3.3 非協調元的分片試驗

為了檢驗采用非協調元的任意網格劃分時是否能達到連續性的要求，需要對其進行分片試驗。若能夠通過分片試驗，則說明解的收斂性能夠得到保證。

(19)

4 有限元實例分析

以單片變截面鋼板彈簧為例，對其進行有限元模型的建立與分析，并將分析結果與其他分析方法比較，從而得出文中分析方法的可行之處。

單片簧的參數包括板簧的展開長度L=1 250 mm、根部截面長度L1=120 mm、端部截面長度L3=220 mm、根部厚度h2=17 mm、端部厚度h1=9 mm，板簧的寬度w=60 mm，板簧自由半徑值R=1 400 mm，以及作用力F=4 500 N、泊松比μ=0.3，彈性模量E=206 000 MPa。以上尺寸參數決定了板簧的截面形狀，由于板簧的對稱性，下面以其半長作為研究對象，并根據其結構特點將它等效為根部固定端部加載作用力進行有限元分析，如圖3所示。

圖3 單片變截面鋼板彈簧半展長

(1)板簧結構分析及節點編排

依據變截面簧的結構特點，對六面體進行有限元網格的劃分、坐標系的確定及網格節點的編排等，見圖4。

圖4 變截面簧網格單元及節點編號示意圖

(2)網格單元定義

針對比較規則的鋼板彈簧，只需要八節點的六面體就能很好地描述其三維實體模型，故這里將對空間六面體單元進行網格單元的定義。網格單元如圖4所示。

(3)賦予材料屬性

彈簧的材料屬性為各向同性材料，當材料有彈性對稱軸時，彈性矩陣中的獨立變量數目將減少，只有彈性模量和泊松比，用于集成彈性矩陣。

(4)添加邊界約束

所謂邊界約束，就是在某個方向上限制物體的位移或者受力。對于鋼板彈簧，根據其安裝位置及結構原理，取半邊作為分析對象，可看作是簡支梁，約束根部節點3個方向的位移，使根部節點1、2、7、8處的X、Y、Z方向位移均為0。

(5)集中力加載

集中力的加載位置為端部沿厚度方向的節點，將載荷力及作用方向按照節點數平均分配到各節點處，比如圖5所示的5、11節點處的y方向。

圖5 文中與ABAQUS有限元網格模型

(6)有限元模型的求解

利用文獻[2]所示的高斯積分法或者Iron積分方法，并基于MATLAB編寫有限元模型的計算程序，根據結構參數建立三維六面體網格的有限元模型，最終計算得出有限元模型的剛度以及節點應力值，如圖6、7所示。

由表1、2可見：采用協調單元計算結果與ABAQUS計算結果存在較大誤差，而采用非協調單元計算結果與ABAQUS計算結果吻合，驗證了非協調單元的計算精度。

圖6 文中有限元分析結果

圖7 ABAQUS分析結果

載荷/N端部位移/mm最大主應力/MPaABAQUS225069.08599.20協調單元[3]225060.16499.36非協調單元(文中)225068.53605.29

表2 與ABAQUS結果的誤差分析 %

5 結論

針對空間六面體協調單元進行推導，并通過分析提出該有限元模型的解剛度值偏大的問題。通過引入了非協調項，使得三維八節點單元中的插值函數二次項完全了，這種三維八節點的非協調元在計算中可以達到與三維二十節點協調元同級的計算精度，但是前者的節點數僅僅是后者的2/5，從而使得在有限元的計算分析中，占用最大內存的平衡方程求解運算效率大大提高，而其計算精度也是非常顯著的。

【1】王勖成.有限單元法[M].北京：清華大學出版社,2003(1):209-216.

【2】王勖成.有限單元法基本原理和數值方法[M].北京：清華大學出版社,1997.

【3】徐榮嬌.結構分析的有限元法與MATLAB程序設計[D].杭州：浙江大學建筑工程學院,2005.

【4】焦兆平.非協調四節點平面等參位移元新列式方法[J].計算結構力學及其應用,1996,13(2):147-156.

【5】龍馭球,辛克貴.廣義協調元[J].清華大學土木工程學報,1987,20(1):1-14.

【6】鹿曉陽,劉玉文,許煥然，等.Wilson非協調元的研究與改進[J].力學學報,1989,21(3):379-384.

【7】陳耀明.汽車懸架論文集[M].蘇州：蘇州大學出版社,2009：32-35.

Study on Finite Element Analysis Method of Incompatible Hexahedral

WEI Jianping1,WEI Hongfa1,WEI Zhian2, HUANG Zhenzhi1

(1.KH Automotive Technologies (Liuzhou) Co.,Ltd., Liuzhou Guangxi 545006,China;2.Liuzhou Wuling Liuji Power Co.,Ltd., Liuzhou Guangxi 545005,China)

Based on the plane iso-parametric incompatible element, hexahedral finite element model expression was shown. Because the stiffness value of complete-compatible element was more large than normal, finally small deformation was lead to. For this purpose, a hexahedral incompatible finite element analysis method was established which could satisfy the fragmentation test condition. And then, low efficiency and big amount of computations for the element of iso-parametric-compatible was solved effectively.Finally, computational accuracy of the incompatible element was verified through example.

Iso-parametric element; Incompatible element; Finite element; Hexahedral mesh

2015-10-26

韋建平(1983—)，男，碩士研究生，從事車輛動力學與控制研究。E-mail:weihf15@163.com。

U464.11

A

1674-1986(2016)02-010-05