基于六面體網(wǎng)格的向量式有限元分析及應(yīng)用

， ，

(1.浙江省建筑設(shè)計(jì)研究院，杭州 310006；2.浙江大學(xué) 空間結(jié)構(gòu)研究中心，杭州 310058)

1 引 言

實(shí)體結(jié)構(gòu)是最為通用的結(jié)構(gòu)形式，具有廣泛的工程應(yīng)用。實(shí)體單元與其他單元通過(guò)在表面連接，以模擬任意外形結(jié)構(gòu)在外界作用下的響應(yīng)。根據(jù)單元形式的不同，實(shí)體結(jié)構(gòu)單元主要有四面體和六面體兩類(lèi)，前者適用性強(qiáng)，而后者計(jì)算精度較高[1,2]。Pawlak等[3]推導(dǎo)了含有旋轉(zhuǎn)自由度的四面體單元形式，該單元對(duì)于具有反對(duì)稱(chēng)應(yīng)力和旋轉(zhuǎn)問(wèn)題的結(jié)構(gòu)計(jì)算具有較高精度。Tayloy等[4]提出位移不協(xié)調(diào)的非協(xié)調(diào)六面體單元形式，并通過(guò)等參變換實(shí)現(xiàn)不規(guī)則網(wǎng)格的數(shù)值模擬。Cao等[5]基于Hu-Washizu變分原理提出多變量形式的六面體單元，提高了單元抗網(wǎng)格畸變能力。李宏光等[6]引入體積坐標(biāo)方法，利用其與整體坐標(biāo)的線性關(guān)系，構(gòu)造出對(duì)網(wǎng)格畸變不敏感的六面體單元，以應(yīng)用于復(fù)雜結(jié)構(gòu)形狀。對(duì)于大位移大轉(zhuǎn)動(dòng)問(wèn)題，六面體單元通過(guò)較少的單元網(wǎng)格劃分即可獲得較為精確的結(jié)果。

向量式有限元VFIFE(vector form intrinsic finite element)[7-9]是基于質(zhì)點(diǎn)運(yùn)動(dòng)和向量分析的新型數(shù)值方法。該法基于牛頓經(jīng)典運(yùn)動(dòng)理論，以質(zhì)點(diǎn)系的運(yùn)動(dòng)來(lái)描述結(jié)構(gòu)體系的行為，各質(zhì)點(diǎn)間的運(yùn)動(dòng)相互獨(dú)立并通過(guò)逐步計(jì)算求得整體結(jié)構(gòu)響應(yīng)，而且求解時(shí)不存在迭代和剛度矩陣問(wèn)題。在變形坐標(biāo)系下，剛體位移可能造成的誤差采用逆向運(yùn)動(dòng)方式來(lái)規(guī)避。故該方法在結(jié)構(gòu)大位移大轉(zhuǎn)動(dòng)、機(jī)構(gòu)運(yùn)動(dòng)及不連續(xù)行為的數(shù)值分析領(lǐng)域具有較好的應(yīng)用[10-14]。

本文給出了八節(jié)點(diǎn)六面體等參實(shí)體單元的向量式基本理論，采用參考面的逆向運(yùn)動(dòng)求解節(jié)點(diǎn)純變形，通過(guò)單元形函的虛功方程計(jì)算節(jié)點(diǎn)內(nèi)力，針對(duì)坐標(biāo)模式和內(nèi)力積分模式等特殊問(wèn)題提出了有效的處理方案。編制數(shù)值計(jì)算程序，并通過(guò)結(jié)構(gòu)算例驗(yàn)證理論和程序的正確性和實(shí)用性。

2 六面體網(wǎng)格的向量式理論

向量式有限元是由質(zhì)點(diǎn)和單元組成，結(jié)構(gòu)的質(zhì)量?jī)H分布于質(zhì)點(diǎn)，本文單元為八節(jié)點(diǎn)六面體等參實(shí)體單元。求解步驟為質(zhì)點(diǎn)全位移計(jì)算、單元節(jié)點(diǎn)純變形計(jì)算和單元節(jié)點(diǎn)內(nèi)力計(jì)算，通過(guò)循環(huán)來(lái)達(dá)到逐點(diǎn)逐步求解。分析流程如圖1所示。

圖1 分析流程圖

Fig.1 Analysis flow chart

2.1 質(zhì)點(diǎn)全位移

質(zhì)點(diǎn)的運(yùn)動(dòng)滿(mǎn)足微分方程：

(1)

式中m為質(zhì)量，α為阻尼參數(shù)，x為位移向量，F(xiàn)=fext+fint為合力向量，上標(biāo)ext和int分別表示外力和內(nèi)力。

時(shí)間步內(nèi)的質(zhì)點(diǎn)全位移采用顯式中央差分方法數(shù)值求解運(yùn)動(dòng)方程獲得，為質(zhì)點(diǎn)新坐標(biāo)與初始坐標(biāo)的差值。

2.2 單元節(jié)點(diǎn)純變形

(1) 坐標(biāo)模式的轉(zhuǎn)換

通過(guò)參考平面的逆向運(yùn)動(dòng)，在質(zhì)點(diǎn)全位移中去除剛體位移，以獲得節(jié)點(diǎn)純變形。節(jié)點(diǎn)實(shí)際坐標(biāo)所組成六面體的每個(gè)面均是空間曲面四邊形，采用投影方式將每個(gè)面均投影為平面四邊形，組合后為投影六面體(即坐標(biāo)模式I)，詳見(jiàn)第3.1節(jié)。

六面體單元面i(i=1～6)的節(jié)點(diǎn)j(j=1,2,3,4)的實(shí)際坐標(biāo)向量記為xj，則平均法向量ni為

(2)

ni x(x-xC i)+ni y(y-yC i)+ni z(z-zC i)=0

(3)

(4)

圖2 六面體單元的投影位置

Fig.2 Projection position of the hexahedral element

(2) 剛體運(yùn)動(dòng)的估算

獲得投影六面體后，采用參考平面估算剛體運(yùn)動(dòng)(平移和轉(zhuǎn)動(dòng))，取投影六面體單元的abcd面為參考平面,如圖3所示。

參考平面abcd在t0～t時(shí)段的空間運(yùn)動(dòng)包括平移、轉(zhuǎn)動(dòng)和純變形[15]。

平移取參考面節(jié)點(diǎn)a的全位移ua，則消去平移后的質(zhì)點(diǎn)位移為

(5)

利用參考面估算質(zhì)點(diǎn)的逆向面外轉(zhuǎn)動(dòng)位移

(6)

質(zhì)點(diǎn)i(i=a～h)的逆向面內(nèi)轉(zhuǎn)動(dòng)位移為

(7)

質(zhì)點(diǎn)純變形是通過(guò)在質(zhì)點(diǎn)全位移中消去平移和轉(zhuǎn)動(dòng)(面外和面內(nèi))后獲得，

(8)

2.3 單元節(jié)點(diǎn)內(nèi)力

(1)傳統(tǒng)形函的引入

圖3 投影六面體單元及其參考平面

Fig.3 Projective hexahedral element and its reference plane

將整體坐標(biāo)系轉(zhuǎn)換至變形坐標(biāo)系，即

(9)

節(jié)點(diǎn)的純變形向量為

(10)

節(jié)點(diǎn)的純變形組合向量為

(11)

(12)

(13)

采用雅可比矩陣J將變形坐標(biāo)系轉(zhuǎn)換至局部坐標(biāo)系，參考文獻(xiàn)[16]。

(14)

(15)

(2) 虛功方程的求解

單元變形滿(mǎn)足虛功方程

(16)

式中fi為整體坐標(biāo)系下的節(jié)點(diǎn)內(nèi)力。

由應(yīng)力應(yīng)變描述的變形虛功為

(17)

單元節(jié)點(diǎn)的內(nèi)力向量為

(18)

(19)

將虛擬狀態(tài)時(shí)的節(jié)點(diǎn)內(nèi)力分量轉(zhuǎn)換至實(shí)際狀態(tài)，并經(jīng)正向運(yùn)動(dòng)，求得實(shí)際位置時(shí)的節(jié)點(diǎn)內(nèi)力：

(20)

單元作用于質(zhì)點(diǎn)的內(nèi)力finti通過(guò)節(jié)點(diǎn)內(nèi)力fi反向作用并集成獲得。

3 關(guān)鍵問(wèn)題的處理

3.1 坐標(biāo)模式

(1) 節(jié)點(diǎn)坐標(biāo)模式

單元節(jié)點(diǎn)坐標(biāo)(坐標(biāo)模式I)用于求解時(shí)間步內(nèi)的節(jié)點(diǎn)純變形和內(nèi)力，因而單元各面上節(jié)點(diǎn)均需處于同一平面。

(2) 質(zhì)點(diǎn)坐標(biāo)模式

3.2 內(nèi)力積分

六面體單元內(nèi)部的應(yīng)力應(yīng)變隨位置變化，求解虛功方程時(shí)需采用數(shù)值積分計(jì)算。采用三維高斯積分方案[1]時(shí)，單積分點(diǎn)和8積分點(diǎn)的相關(guān)系數(shù)列入表1，本文采用后者求解式(18)，可獲得較好計(jì)算結(jié)果。

表1 數(shù)值積分方案

Tab.1 Numerical integration programmes

積分點(diǎn)數(shù)焊iηiζiHijm10008.08±1/3±1/3±1/31.0,1.0,1.0,1.01.0,1.0,1.0,1.0

(21)

(22)

式中k為每個(gè)坐標(biāo)方向的積分點(diǎn)數(shù)。

4 工程結(jié)構(gòu)算例分析

基于向量式有限元六面體等參元理論，在Matlab平臺(tái)上編制分析計(jì)算程序，實(shí)現(xiàn)向量式有限元實(shí)體結(jié)構(gòu)的計(jì)算分析。

4.1 懸臂梁結(jié)構(gòu)靜動(dòng)力分析

(1) 靜力分析

懸臂梁頂面施加豎向均布靜力作用q=2.5×104kN/m2，采用斜坡-平臺(tái)方式緩慢加載，并考慮阻尼效應(yīng)來(lái)獲得收斂值。細(xì)長(zhǎng)梁自由端的撓度理論解為

w(x)=(x2-4Lx+6L2)(qH)x2/(24EI)

圖5為懸臂梁自由端頂部節(jié)點(diǎn)豎向位移的計(jì)算收斂過(guò)程，圖6為沿懸臂梁頂部中線各節(jié)點(diǎn)的豎向位移和von Mises應(yīng)力結(jié)果?？梢钥闯觯Q向位移和von Mises應(yīng)力結(jié)果與細(xì)長(zhǎng)梁理論解及ABAQUS結(jié)果均基本一致。本文和ABAQUS的自由端節(jié)點(diǎn)豎向位移分別為0.0115 m和0.0119 m，相對(duì)理論解0.0117 m的誤差分別為-1.71%和1.71%。本文和ABAQUS的自由端節(jié)點(diǎn)von Mises應(yīng)力分別為38.06 MPa和36.4 MPa，差異約為4.36%，懸臂端節(jié)點(diǎn)由于應(yīng)力集中，結(jié)果差異相對(duì)較大，約為9.8%?？梢?jiàn)本文方法在實(shí)體工程結(jié)構(gòu)靜力分析中具有較好準(zhǔn)確性。

(2) 動(dòng)力分析

在懸臂梁頂面施加瞬時(shí)動(dòng)力均布荷載q=2.0×104kN/m2，不考慮阻尼效應(yīng)。取t=0.01 s進(jìn)行分析，圖7和圖8分別為自由端頂部節(jié)點(diǎn)豎向位移和von Mises應(yīng)力隨時(shí)間的變化情況。

結(jié)果表明，本文位移和應(yīng)力結(jié)果均與ABAQUS結(jié)果吻合較好，驗(yàn)證了本文方法在實(shí)體工程結(jié)構(gòu)動(dòng)力分析中的準(zhǔn)確性。圖8中兩者應(yīng)力結(jié)果的總體變化趨勢(shì)基本一致，由于端部應(yīng)力集中，局部振蕩波形有所差異。

圖4 懸臂梁模型

Fig.4 Cantilever-beam model

圖5 自由端豎向位移收斂過(guò)程

Fig.5 Convergence process of free end vertical displacement

圖6 頂部中線位置結(jié)果

Fig.6 Results of the top center-line

4.2 環(huán)形套筒大變形分析

取t=0.02 s進(jìn)行分析。圖10為環(huán)形套筒左側(cè)中心節(jié)點(diǎn)的水平位移隨時(shí)間的變化情況?？梢钥闯觯疚姆椒色@得套筒大變形過(guò)程中的位移變化情況，且與ABAQUS結(jié)果基本一致。在同時(shí)采用顯式動(dòng)力分析時(shí)，本文方法由于不存在剛度矩陣奇異和迭代不收斂問(wèn)題，相比ABAQUS具有更高的計(jì)算效率。

典型的環(huán)形套筒變形和von Mises應(yīng)力云圖情況如圖11所示?？梢钥闯?，套筒起初表現(xiàn)為壓扁變形(t=0.004 s)，t=0.008 s時(shí)，壓扁變形達(dá)到最大值；接著套筒開(kāi)始回彈，在t=0.012 s時(shí)回彈至最小變形；然后，套筒壓扁變形再次變大(t=0.015 s)；依此循環(huán)往復(fù)。本文所得結(jié)構(gòu)變形和應(yīng)力分布情況均與ABAQUS計(jì)算結(jié)果吻合。

圖7 自由端豎向位移-時(shí)間曲線

Fig.7 Vertical displacement-time curves of the free end

圖8 自由端von Mises應(yīng)力-時(shí)間曲線

Fig.8 Von Mises stress-time curves of the free end

圖9 環(huán)形套筒模型

Fig.9 Annular sleeve model

圖10 左側(cè)中心節(jié)點(diǎn)的位移-時(shí)間曲線

Fig.10 Displacement-time curves of the left central node

圖11 環(huán)形套筒壓扁變形圖

Fig.11 Flattened deformation diagrams of annular sleeve

5 結(jié) 論

基于向量式有限元，給出了八節(jié)點(diǎn)六面體等參實(shí)體單元的基本公式，通過(guò)投影方式將空間曲面六面體轉(zhuǎn)換為投影六面體，采用參考面的逆向運(yùn)動(dòng)求解節(jié)點(diǎn)純變形，通過(guò)單元形函的虛功方程計(jì)算節(jié)點(diǎn)內(nèi)力。針對(duì)坐標(biāo)模式和內(nèi)力積分模式等關(guān)鍵問(wèn)題提出了有效的處理方案。

編寫(xiě)了六面體實(shí)體單元的數(shù)值計(jì)算程序，并針對(duì)工程結(jié)構(gòu)的懸臂梁靜動(dòng)力和環(huán)形套筒大變形算例進(jìn)行分析，驗(yàn)證了本文理論和程序的正確性和實(shí)用性。

向量式有限元適用于分析工程結(jié)構(gòu)大變形大位移行為的動(dòng)力分析，求解時(shí)不存在迭代和剛度矩陣問(wèn)題，在結(jié)構(gòu)動(dòng)力分析領(lǐng)域具有更高的效率和準(zhǔn)確性。