非協(xié)調(diào)六面體有限元分析方法研究

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

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

?

非協(xié)調(diào)六面體有限元分析方法研究

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

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

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

等參元；非協(xié)調(diào)元；有限元；六面體網(wǎng)格

0 引言

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

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

1 有限元力學(xué)基礎(chǔ)

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

(1)

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

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

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

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

(2)

nσ=ps(在邊界V上)

2 協(xié)調(diào)六面體等參有限元

2.1 六面體單元概述

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

圖1 四邊形等參數(shù)單元

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

圖2 八節(jié)點(diǎn)六面體等參單元

2.2 六面體單元特性

根據(jù)等參數(shù)單元的含義，它的位移和坐標(biāo)都采用相同的形函數(shù)表示，即

位移函數(shù)：

(3)

坐標(biāo)函數(shù)：

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

(4)

式中：ξi，ηi和ζi是節(jié)點(diǎn)i的局部坐標(biāo)，分別位于六面體的8個角處，其對應(yīng)坐標(biāo)值分別為1或者-1，即：

(5)

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

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

(6)

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

(7)

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

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

(8)

式中：B稱為單元應(yīng)變系數(shù)，從而單元應(yīng)力的計算公式為：

(9)

式中：

(10)

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

(11)

綜上所述，并根據(jù)單元節(jié)點(diǎn)的剛度公式，三維等參單元的剛度矩陣表達(dá)式如下：

(12)

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

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

(13)

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

3 非協(xié)調(diào)元的凝聚

3.1 Wilson非協(xié)調(diào)元

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

而附加項的位移在三維線性單元的8個節(jié)點(diǎn)上都取零值時，它對節(jié)點(diǎn)位移沒有任何影響，而只對單元內(nèi)部的位移起了調(diào)整作用。待定系數(shù)α1、α2和α3與單元邊界上的節(jié)點(diǎn)無關(guān)。

包含附加的無節(jié)點(diǎn)位移項的單元位移插值表示如下：

(14)

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

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

因此，為了改善三維線性單元的性質(zhì)，解決完備協(xié)調(diào)元中剛度系數(shù)值過大，提高其解精度，需將單元的非協(xié)調(diào)位移插值函數(shù)附加到內(nèi)部無節(jié)點(diǎn)的位移項。

3.2 非協(xié)調(diào)元剛度的凝聚

根據(jù)上文所述，結(jié)合幾何關(guān)系和位移函數(shù)，加入非協(xié)調(diào)部分后，單元的應(yīng)變矩陣公式變?yōu)槿缦率阶樱?/p>

(15)

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

(16)

其中：

(17)

式中：Pe為外載荷列陣。

對上述關(guān)系式進(jìn)行展開求解，則單元內(nèi)部位移、剛度矩陣以及載荷列陣之間的關(guān)系式

Keδe=Pe

(18)

式中：

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

3.3 非協(xié)調(diào)元的分片試驗

為了檢驗采用非協(xié)調(diào)元的任意網(wǎng)格劃分時是否能達(dá)到連續(xù)性的要求，需要對其進(jìn)行分片試驗。若能夠通過分片試驗，則說明解的收斂性能夠得到保證。

(19)

4 有限元實例分析

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

單片簧的參數(shù)包括板簧的展開長度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。以上尺寸參數(shù)決定了板簧的截面形狀，由于板簧的對稱性，下面以其半長作為研究對象，并根據(jù)其結(jié)構(gòu)特點(diǎn)將它等效為根部固定端部加載作用力進(jìn)行有限元分析，如圖3所示。

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

(1)板簧結(jié)構(gòu)分析及節(jié)點(diǎn)編排

依據(jù)變截面簧的結(jié)構(gòu)特點(diǎn)，對六面體進(jìn)行有限元網(wǎng)格的劃分、坐標(biāo)系的確定及網(wǎng)格節(jié)點(diǎn)的編排等，見圖4。

圖4 變截面簧網(wǎng)格單元及節(jié)點(diǎn)編號示意圖

(2)網(wǎng)格單元定義

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

(3)賦予材料屬性

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

(4)添加邊界約束

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

(5)集中力加載

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

圖5 文中與ABAQUS有限元網(wǎng)格模型

(6)有限元模型的求解

利用文獻(xiàn)[2]所示的高斯積分法或者Iron積分方法，并基于MATLAB編寫有限元模型的計算程序，根據(jù)結(jié)構(gòu)參數(shù)建立三維六面體網(wǎng)格的有限元模型，最終計算得出有限元模型的剛度以及節(jié)點(diǎn)應(yīng)力值，如圖6、7所示。

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

圖6 文中有限元分析結(jié)果

圖7 ABAQUS分析結(jié)果

載荷/N端部位移/mm最大主應(yīng)力/MPaABAQUS225069.08599.20協(xié)調(diào)單元[3]225060.16499.36非協(xié)調(diào)單元(文中)225068.53605.29

表2 與ABAQUS結(jié)果的誤差分析 %

5 結(jié)論

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

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

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

【3】徐榮嬌.結(jié)構(gòu)分析的有限元法與MATLAB程序設(shè)計[D].杭州：浙江大學(xué)建筑工程學(xué)院,2005.

【4】焦兆平.非協(xié)調(diào)四節(jié)點(diǎn)平面等參位移元新列式方法[J].計算結(jié)構(gòu)力學(xué)及其應(yīng)用,1996,13(2):147-156.

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

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

【7】陳耀明.汽車懸架論文集[M].蘇州：蘇州大學(xué)出版社,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—)，男，碩士研究生，從事車輛動力學(xué)與控制研究。E-mail:weihf15@163.com。

U464.11

A

1674-1986(2016)02-010-05