, ,
(1.浙江省建筑設(shè)計(jì)研究院,杭州 310006;2.浙江大學(xué) 空間結(jié)構(gòu)研究中心,杭州 310058)
實(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í)用性。
向量式有限元是由質(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
質(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)的差值。
(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)
(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反向作用并集成獲得。
(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)模式
六面體單元內(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ù)。
基于向量式有限元六面體等參元理論,在Matlab平臺(tái)上編制分析計(jì)算程序,實(shí)現(xiàn)向量式有限元實(shí)體結(jié)構(gòu)的計(jì)算分析。
(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
取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
基于向量式有限元,給出了八節(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)確性。