， ，

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

1 引 言

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

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

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

2 六面體網格的向量式理論

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

圖1 分析流程圖

Fig.1 Analysis flow chart

2.1 質點全位移

質點的運動滿足微分方程：

(1)

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

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

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

(1) 坐標模式的轉換

通過參考平面的逆向運動，在質點全位移中去除剛體位移，以獲得節(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) 剛體運動的估算

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

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

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

(5)

利用參考面估算質點的逆向面外轉動位移

(6)

質點i(i=a～h)的逆向面內轉動位移為

(7)

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

(8)

2.3 單元節(jié)點內力

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

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

Fig.3 Projective hexahedral element and its reference plane

將整體坐標系轉換至變形坐標系，即

(9)

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

(10)

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

(11)

(12)

(13)

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

(14)

(15)

(2) 虛功方程的求解

單元變形滿足虛功方程

(16)

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

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

(17)

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

(18)

(19)

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

(20)

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

3 關鍵問題的處理

3.1 坐標模式

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

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

(2) 質點坐標模式

3.2 內力積分

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

表1 數值積分方案

Tab.1 Numerical integration programmes

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

4 工程結構算例分析

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

4.1 懸臂梁結構靜動力分析

(1) 靜力分析

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

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

圖5為懸臂梁自由端頂部節(jié)點豎向位移的計算收斂過程，圖6為沿懸臂梁頂部中線各節(jié)點的豎向位移和von Mises應力結果?？梢钥闯觯Q向位移和von Mises應力結果與細長梁理論解及ABAQUS結果均基本一致。本文和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é)點由于應力集中，結果差異相對較大，約為9.8%?？梢姳疚姆椒ㄔ趯嶓w工程結構靜力分析中具有較好準確性。

(2) 動力分析

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

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

圖4 懸臂梁模型

Fig.4 Cantilever-beam model

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

Fig.5 Convergence process of free end vertical displacement

圖6 頂部中線位置結果

Fig.6 Results of the top center-line

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

取t=0.02 s進行分析。圖10為環(huán)形套筒左側中心節(jié)點的水平位移隨時間的變化情況?？梢钥闯觯疚姆椒色@得套筒大變形過程中的位移變化情況，且與ABAQUS結果基本一致。在同時采用顯式動力分析時，本文方法由于不存在剛度矩陣奇異和迭代不收斂問題，相比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)往復。本文所得結構變形和應力分布情況均與ABAQUS計算結果吻合。

圖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 左側中心節(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é)點純變形，通過單元形函的虛功方程計算節(jié)點內力。針對坐標模式和內力積分模式等關鍵問題提出了有效的處理方案。

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

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