， ，

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

1 引 言

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

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

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

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

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

圖1 分析流程圖

Fig.1 Analysis flow chart

2.1 質(zhì)點全位移

質(zhì)點的運動滿足微分方程：

(1)

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

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

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

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

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

六面體單元面i(i=1～6)的節(jié)點j(j=1,2,3,4)的實際坐標向量記為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) 剛體運動的估算

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

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

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

(5)

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

(6)

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

(7)

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

(8)

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

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

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

Fig.3 Projective hexahedral element and its reference plane

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

(9)

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

(10)

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

(11)

(12)

(13)

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

(14)

(15)

(2) 虛功方程的求解

單元變形滿足虛功方程

(16)

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

由應力應變描述的變形虛功為

(17)

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

(18)

(19)

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

(20)

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

3 關鍵問題的處理

3.1 坐標模式

(1) 節(jié)點坐標模式

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

(2) 質(zhì)點坐標模式

3.2 內(nèi)力積分

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

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

Tab.1 Numerical integration programmes

積分點數(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為每個坐標方向的積分點數(shù)。

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

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

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

(1) 靜力分析

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

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

圖5為懸臂梁自由端頂部節(jié)點豎向位移的計算收斂過程，圖6為沿懸臂梁頂部中線各節(jié)點的豎向位移和von Mises應力結(jié)果?？梢钥闯?，豎向位移和von Mises應力結(jié)果與細長梁理論解及ABAQUS結(jié)果均基本一致。本文和ABAQUS的自由端節(jié)點豎向位移分別為0.0115 m和0.0119 m，相對理論解0.0117 m的誤差分別為-1.71%和1.71%。本文和ABAQUS的自由端節(jié)點von Mises應力分別為38.06 MPa和36.4 MPa，差異約為4.36%，懸臂端節(jié)點由于應力集中，結(jié)果差異相對較大，約為9.8%?？梢姳疚姆椒ㄔ趯嶓w工程結(jié)構(gòu)靜力分析中具有較好準確性。

(2) 動力分析

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

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

圖4 懸臂梁模型

Fig.4 Cantilever-beam model

圖5 自由端豎向位移收斂過程

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進行分析。圖10為環(huán)形套筒左側(cè)中心節(jié)點的水平位移隨時間的變化情況。可以看出，本文方法可獲得套筒大變形過程中的位移變化情況，且與ABAQUS結(jié)果基本一致。在同時采用顯式動力分析時，本文方法由于不存在剛度矩陣奇異和迭代不收斂問題，相比ABAQUS具有更高的計算效率。

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

圖7 自由端豎向位移-時間曲線

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

圖8 自由端von Mises應力-時間曲線

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

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

Fig.9 Annular sleeve model

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

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

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

Fig.11 Flattened deformation diagrams of annular sleeve

5 結(jié) 論

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

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

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