王東東, 吳振宇, 侯松陽

(廈門大學(xué) 土木工程系 福建省濱海土木工程數(shù)字仿真重點(diǎn)實(shí)驗(yàn)室,廈門 361005)

1 引 言

工程中忽略剪切變形的細(xì)長梁結(jié)構(gòu)常稱作歐拉伯努利梁或歐拉梁[1]。由于歐拉梁問題的控制方程為四階微分方程,其等效積分弱形式包含位移的二階導(dǎo)數(shù),因此要求數(shù)值方法的形函數(shù)至少具有C1連續(xù)性[2-4]。在有限元法中,應(yīng)用最為廣泛的歐拉梁單元是兩節(jié)點(diǎn)三次埃爾米特梁單元[3,4]。該單元兩個(gè)節(jié)點(diǎn)均包含撓度和轉(zhuǎn)角自由度,在單元邊界處滿足C1連續(xù)性要求,而且撓度和轉(zhuǎn)角自由度對應(yīng)的形函數(shù)滿足埃爾米特插值特性,可以便捷施加撓度和轉(zhuǎn)角邊界條件。另一方面,無網(wǎng)格法和等幾何分析采用的形函數(shù)具有任意高階連續(xù)的特性,但這類形函數(shù)通常沒有插值特性[5-9]。對于一維等幾何開放型B樣條基函數(shù),其在兩個(gè)端點(diǎn)滿足插值特性,但是由于不包含轉(zhuǎn)角自由度,不具有轉(zhuǎn)角插值性。因此,無網(wǎng)格法和等幾何分析一般需要引入其他方法來處理邊界條件[6-9]。

在結(jié)構(gòu)動力有限元分析中,質(zhì)量矩陣的構(gòu)造對計(jì)算精度有明顯的影響。常用的質(zhì)量矩陣包括一致質(zhì)量矩陣[3]、集中質(zhì)量矩陣[10-15]和高階質(zhì)量矩陣[3,16-18]。其中,集中質(zhì)量矩陣構(gòu)造簡單,與顯式時(shí)間積分結(jié)合使用可以大幅提高結(jié)構(gòu)動力計(jì)算效率。對于三次埃爾米特梁單元,其集中質(zhì)量矩陣通常不考慮轉(zhuǎn)角慣性矩的影響,而由撓度自由度對應(yīng)的一致質(zhì)量矩陣元素通過行求和或節(jié)點(diǎn)積分構(gòu)造[3]。對于兩端簡支或固定的梁問題,該集中質(zhì)量矩陣的頻率解呈現(xiàn)最優(yōu)的4階精度[3]。但對包含自由端的梁問題,采用該集中質(zhì)量矩陣計(jì)算得到的頻率精度僅為2階,出現(xiàn)明顯的精度掉階問題。此外,目前文獻(xiàn)中討論的埃爾米特梁單元主要是指兩節(jié)點(diǎn)三次單元,還缺乏高次單元的統(tǒng)一構(gòu)造方式。最近,Wu等[18]在埃爾米特再生核無網(wǎng)格法的理論框架內(nèi)[19-21],從理論上證明當(dāng)埃爾米特再生核無網(wǎng)格撓度和轉(zhuǎn)角形函數(shù)的影響域與單元完全重合時(shí),其退化為埃爾米特有限元撓度和轉(zhuǎn)角形函數(shù),從而為任意次埃爾米特單元的撓度和轉(zhuǎn)角形函數(shù)構(gòu)造提供了一種統(tǒng)一方法。根據(jù)該理論,可以方便地構(gòu)造三節(jié)點(diǎn)五次埃爾米特單元。對于五次埃爾米特單元,當(dāng)采用忽略旋轉(zhuǎn)影響的行求和質(zhì)量矩陣進(jìn)行包含自由邊界梁問題振動分析時(shí),頻率精度同樣只有2階,遠(yuǎn)低于理論收斂階次2(p-1)[22],其中p為單元的階次。由此可見,如何構(gòu)造能有效規(guī)避邊界條件影響的集中質(zhì)量矩陣是埃爾米特梁單元的一個(gè)共性問題。

針對上述問題,本文從質(zhì)量矩陣的數(shù)值積分精度出發(fā),詳細(xì)討論了滿足理論積分精度要求的梯度增強(qiáng)節(jié)點(diǎn)積分方案,厘清了集中質(zhì)量矩陣頻率精度掉階的原因:從節(jié)點(diǎn)積分角度而言,質(zhì)量矩陣的數(shù)值積分精度不滿足理論最低積分精度的要求。隨后,根據(jù)不同單元質(zhì)量矩陣的數(shù)值積分精度要求,發(fā)展了三次和五次單元梯度增強(qiáng)節(jié)點(diǎn)積分方案,并以此為基礎(chǔ)構(gòu)造了相應(yīng)的分塊對角形式的單元質(zhì)量矩陣,其組裝得到的整體質(zhì)量矩陣除邊界節(jié)點(diǎn)對應(yīng)的元素外呈現(xiàn)對角形式。其中,三次單元的分塊對角矩陣不但具有與一致質(zhì)量矩陣相同的最優(yōu)收斂特性,而且不依賴于邊界條件;五次單元的分塊對角矩陣的精度同樣不依賴于邊界條件,只是精度呈現(xiàn)次優(yōu)收斂性。為了進(jìn)一步提升頻率計(jì)算精度,通過將一致質(zhì)量矩陣和基于節(jié)點(diǎn)積分質(zhì)量矩陣進(jìn)行組合,建立了三次和五次單元的高階質(zhì)量矩陣,其精度比對應(yīng)的一致質(zhì)量矩陣高2階,同時(shí)不受邊界條件和波數(shù)等因素的影響。

2 三次和五次梁單元與集中質(zhì)量矩陣

2.1 三次和五次埃爾米特梁單元

考慮空間區(qū)域Ω的歐拉梁問題,其自由振動問題的有限元離散方程為

[K-(ωh)2M]d=0

(1)

式中ωh為頻率數(shù)值解,d為相應(yīng)的振動模態(tài),M和K分別為整體質(zhì)量矩陣和整體剛度矩陣,可由相應(yīng)的單元質(zhì)量矩陣和單元剛度矩陣組裝得到

(2)

圖1 三次梁單元的撓度及轉(zhuǎn)角形函數(shù)

圖2 五次梁單元的撓度及轉(zhuǎn)角形函數(shù)

(3)

(1) 三次單元

(4)

(5)

(6)

(7)

(2) 五次單元

(8)

(9)

(10)

(11)

2.2 集中質(zhì)量矩陣的精度掉階問題

由圖3和圖4可以看出,采用一致質(zhì)量矩陣計(jì)算的頻率結(jié)果均能達(dá)到最優(yōu)收斂率,即三次單元為4階,五次單元為8階。對于簡支梁,三次單元集中質(zhì)量矩陣同樣能夠達(dá)到4階最優(yōu)收斂,而五次單元集中質(zhì)量矩陣的頻率則是6階收斂,比最優(yōu)收斂率低2階。但對于兩端自由梁,無論三次單元還是五次單元,集中質(zhì)量矩陣的頻率收斂率都驟降至2階,遠(yuǎn)低于最優(yōu)收斂率,呈現(xiàn)明顯的頻率精度掉階問題。本文從質(zhì)量矩陣的節(jié)點(diǎn)積分方案入手,揭示頻率精度掉階的內(nèi)在機(jī)理。

圖3 簡支梁的基頻收斂結(jié)果對比

圖4 兩端自由梁的基頻收斂結(jié)果對比

3 三次梁單元分塊對角與高階質(zhì)量矩陣

3.1 三次埃爾米特梁單元的分塊對角質(zhì)量矩陣

傳統(tǒng)的節(jié)點(diǎn)積分方案將單元的節(jié)點(diǎn)作為積分點(diǎn),相應(yīng)質(zhì)量矩陣Me的計(jì)算式為

(12)

(13)

式中m=0,1,…,nen-1。

此外,注意到埃爾米特形函數(shù)及其導(dǎo)數(shù)滿足如下的插值性質(zhì)[18]

(14)

(15)

式中δij為克羅內(nèi)克符號。傳統(tǒng)的節(jié)點(diǎn)積分法一般僅利用式(14)的撓度形函數(shù)插值特性。根據(jù)式(12,13),三次單元節(jié)點(diǎn)積分方案為常見的矩形法則

{ξi}={-1,1}, {wi}={1,1}

(16)

根據(jù)式(16)得到的Me即為式(7)的三次單元集中質(zhì)量矩陣。值得指出的是,集中質(zhì)量矩陣若要達(dá)到與一致質(zhì)量矩陣相同的最優(yōu)收斂特性,計(jì)算集中質(zhì)量矩陣所用節(jié)點(diǎn)積分方案應(yīng)至少具有2(p-2)階精度[3,22]。由于式(16)的三次單元節(jié)點(diǎn)積分方案僅有1階精度,低于2階理論精度,故相應(yīng)的集中質(zhì)量矩陣給出的頻率達(dá)不到最優(yōu)收斂率。

為了提高節(jié)點(diǎn)積分方案的精度,可利用埃爾米特形函數(shù)C1連續(xù)的優(yōu)點(diǎn),構(gòu)造如下的梯度增強(qiáng)節(jié)點(diǎn)積分方案

(17)

(18)

式中m=0,1,…,2nen-1。

根據(jù)式(17,18),可得具有3階精度的兩節(jié)點(diǎn)三次單元梯度增強(qiáng)節(jié)點(diǎn)積分方案,即

(19)

將式(19)代入式(17),有

(20)

如前所述,該質(zhì)量矩陣滿足數(shù)值積分精度要求,因此對應(yīng)的頻率解具有4階最優(yōu)收斂率。值得注意的是,如式(14,15)所示,埃爾米特?fù)隙刃魏瘮?shù)和轉(zhuǎn)角形函數(shù)的導(dǎo)數(shù)在單元節(jié)點(diǎn)上均具有插值性,因此采用梯度增強(qiáng)節(jié)點(diǎn)積分法,根據(jù)式(17,18)計(jì)算得到的單元質(zhì)量矩陣具有分塊對角形式,本文稱之為分塊對角質(zhì)量矩陣。

3.2 三次埃爾米特梁單元的高階質(zhì)量矩陣

(21)

為了衡量頻率計(jì)算精度,引入如下相對誤差

e=ωh/ω-1,ωh=(1+e)ω

(22)

考慮如圖5所示的典型三次單元均勻有限元網(wǎng)格,引入如下的節(jié)點(diǎn)位移簡諧波假定

(23)

(24)

式中

A11=[-54(ωh)2h4β-5040c2]cos(kh)+54(ωh)2h4(β-35/9)+5040c2

(25)

(26)

A22=3h[(ωh)2h4β+280c2]cos(kh)-

4h[(ωh)2h4β-420c2]

(27)

圖5 三次埃爾米特單元節(jié)點(diǎn)位移的簡諧波表示

由式(24)的非零解條件,可知det[Aij]=0,同時(shí)利用式(22),并忽略誤差e高階項(xiàng),有

e≈-a0/a1

(28)

式中

1440(β+35/2)k4h4]cos(kh)+302400+(β2k8h8+720βk4h4+302400)cos(kh)2-

55β(β-24/11)k8h8+720(β-70)k4h4}

(29)

[90cos(kh)+55β+120]βk4h4+360cos(kh)[βcos(kh)-2β-35]+360(β-70)}

(30)

進(jìn)一步對式(28)的余弦項(xiàng)關(guān)于h進(jìn)行泰勒展開,可得Me(β)頻率誤差表達(dá)式為

(31)

顯然,當(dāng)β=1/2時(shí)可獲得6階超收斂結(jié)果。此時(shí)的高階質(zhì)量矩陣以及對應(yīng)的頻率精度分別為

(32)

(33)

(34)

(35)

其中,式(34)的誤差結(jié)果從理論上驗(yàn)證了分塊對角質(zhì)量矩陣得到的頻率具有4階最優(yōu)收斂性的結(jié)論。

4 五次梁單元分塊對角與高階質(zhì)量矩陣

4.1 五次埃爾米特梁單元的分塊對角質(zhì)量矩陣

根據(jù)式(18),經(jīng)計(jì)算可得三節(jié)點(diǎn)五次單元的梯度增強(qiáng)節(jié)點(diǎn)積分方案,即

(36)

該積分方案具有5階精度,對應(yīng)的質(zhì)量矩陣為

(37)

(38)

(39)

(40)

(41)

將式(41)代入式(39),有

(42)

4.2 五次埃爾米特梁單元的高階質(zhì)量矩陣

(43)

對其進(jìn)行頻率精度分析,可得

k10h10+O(h12)

(44)

顯然,當(dāng)γ=13/21時(shí),頻率精度為10階,具有超收斂特性,相應(yīng)的高階質(zhì)量矩陣以及頻率誤差為

(45)

(46)

(47)

(48)

5 數(shù)值驗(yàn)證

為了全面驗(yàn)證本文分塊對角質(zhì)量矩陣的魯棒性與高階質(zhì)量矩陣的超收斂特性,采用18,26和38個(gè)三次單元和9,13和19個(gè)五次單元來計(jì)算圖6所示的四種不同邊界條件下梁的振動問題,即兩端固支梁、簡支梁、懸臂梁和兩端自由梁,梁的幾何和材料參數(shù)已在2.2節(jié)給出。由于基頻的計(jì)算精度接近機(jī)器精度,為了更清楚地說明所提方法的收斂率和精度,圖6給出了采用不同質(zhì)量矩陣計(jì)算得到的第七階頻率的收斂結(jié)果對比。結(jié)果表明,傳統(tǒng)的集中質(zhì)量矩陣在含有自由端邊界條件下會出現(xiàn)嚴(yán)重的頻率精度掉階問題。與此同時(shí),基于梯度增強(qiáng)節(jié)點(diǎn)積分法的對角分塊質(zhì)量矩陣具有良好的魯棒性,相應(yīng)三次和五次單元的頻率計(jì)算結(jié)果分別達(dá)到4階最優(yōu)收斂和6階次優(yōu)收斂。此外,由一致質(zhì)量矩陣和基于梯度增強(qiáng)節(jié)點(diǎn)積分方案構(gòu)建的質(zhì)量矩陣進(jìn)行最優(yōu)線性組合得到的高階質(zhì)量矩陣同樣不受邊界條件的影響,對于三次和五次單元,頻率結(jié)果呈現(xiàn)出6階和10階的超收斂特性,分別比相應(yīng)的一致質(zhì)量矩陣頻率精度高出2階。

圖6 不同邊界條件梁問題第七階頻率收斂結(jié)果對比

6 結(jié) 論

本文針對埃爾米特梁單元集中質(zhì)量矩陣受邊界條件影響的頻率精度掉階問題,從質(zhì)量矩陣構(gòu)造過程中的數(shù)值積分精度進(jìn)行分析,闡明了精度掉階的原因在于單元集中質(zhì)量矩陣的構(gòu)造積分方法不滿足最優(yōu)收斂性對應(yīng)的積分精度要求。在此基礎(chǔ)上,利用埃爾米特梁單元撓度形函數(shù)和轉(zhuǎn)角形函數(shù)導(dǎo)數(shù)均具有節(jié)點(diǎn)插值特性,發(fā)展了考慮梯度項(xiàng)的梯度增強(qiáng)節(jié)點(diǎn)積分方案,進(jìn)而建立了具有分塊對角形式的單元質(zhì)量矩陣。具體而言,三次和五次單元的分塊對角質(zhì)量矩陣的頻率解分別具有穩(wěn)定的4階最優(yōu)收斂和6階次優(yōu)收斂性質(zhì)。需要指出的是,對于均布單元,由單元分塊質(zhì)量矩陣組裝得到的整體質(zhì)量矩陣除邊界節(jié)點(diǎn)相關(guān)元素外仍然保持對角形式。為了進(jìn)一步提高頻率計(jì)算精度,通過將一致質(zhì)量矩陣和與其具有同階精度的梯度增強(qiáng)節(jié)點(diǎn)積分質(zhì)量矩陣進(jìn)行優(yōu)化組合,建立了不受邊界條件、波數(shù)等因素影響的高階質(zhì)量矩陣,相應(yīng)頻率精度比一致質(zhì)量矩陣提升了2階,具有超收斂特性。不同邊界條件下梁的振動分析結(jié)果顯示,埃爾米特梁單元的分塊對角質(zhì)量矩陣和高階質(zhì)量矩陣的數(shù)值頻率均能達(dá)到理論精度,并具有良好的魯棒性。