空間剛架穩(wěn)定性分析的泡函數(shù)有限元法

陳 梁，陳振中(沈陽(yáng)航空航天大學(xué) 航空航天工程學(xué)部(院)，沈陽(yáng) 110136)

空間剛架穩(wěn)定性分析的泡函數(shù)有限元法

陳 梁，陳振中
(沈陽(yáng)航空航天大學(xué) 航空航天工程學(xué)部(院)，沈陽(yáng) 110136)

將泡函數(shù)有限元法理論用于空間剛架的穩(wěn)定性分析。為了全面地考慮各種載荷因素對(duì)幾何剛度矩陣的影響，提出了變幾何剛度法和失穩(wěn)判據(jù)。為了解決單元坐標(biāo)系方位的確定這一個(gè)性化情況很突出的問(wèn)題，運(yùn)用三點(diǎn)定位算法編制了規(guī)范的計(jì)算程序。在一桿一單元的情況下，運(yùn)用這三大理論編程計(jì)算出了比較精確的臨界載荷值，說(shuō)明在空間剛架的穩(wěn)定性計(jì)算方面，泡函數(shù)有限元法具有突出的簡(jiǎn)潔性和精確性。

空間剛架結(jié)構(gòu)；穩(wěn)定性分析；泡函數(shù)；變幾何剛度法；三點(diǎn)定位法；有限單元法

在生產(chǎn)實(shí)踐中存在大量的空間剛架結(jié)構(gòu)，如塔吊、汽車吊、輸變電塔、鋼構(gòu)橋等，在遭受橫向擾動(dòng)的情況下這些剛架能否做到不失穩(wěn)將是考驗(yàn)該類產(chǎn)品設(shè)計(jì)成敗的關(guān)鍵。然而這類剛架存在桿件數(shù)量多，桿件截面空間布置的個(gè)性化突出等實(shí)際情況，這就使空間復(fù)雜剛架系統(tǒng)的穩(wěn)定性分析具有很大的挑戰(zhàn)性。王仁輝、高軒能等[1]將泡函數(shù)引入平面剛架的穩(wěn)定性計(jì)算，并獲得了簡(jiǎn)潔而精確的結(jié)果；龍馭球、袁駟等[2]提出極值點(diǎn)失穩(wěn)問(wèn)題的處理方案。在前人研究成果的基礎(chǔ)上，本文將泡函數(shù)有限元法應(yīng)用于空間剛架穩(wěn)定性分析。為了全面地考慮各種載荷因素對(duì)幾何剛度矩陣的影響，提出了變幾何剛度法和失穩(wěn)判據(jù)。為了解決繁瑣的桿件單元坐標(biāo)系的確定問(wèn)題，利用三點(diǎn)定位算法編制了規(guī)范的計(jì)算程序。

1 泡函數(shù)有限元法

在單一載荷(即待求的臨界載荷)作用下，結(jié)構(gòu)的失穩(wěn)形式為分岔點(diǎn)失穩(wěn)。采用一桿一單元，以桿件的一端為原點(diǎn)建立單元局部坐標(biāo)系o-x′y′，x′軸的正向指向桿的另一端，在剛架所在平面內(nèi)y′軸垂直于x′指向桿件的一側(cè)。 設(shè)桿件的長(zhǎng)度為L(zhǎng)，橫截面積為A，彈性模量為E，抗彎軸慣性矩為I。根據(jù)文獻(xiàn)[1]，只考慮彎曲變形影響的泡函數(shù)單元的形函數(shù)矩陣為

N=[N1，N2，N3，N4，N5],N1=(2η+1)(η-1)2,N2=Lη(η-1)2,N3=-η2(2η-3),N4=Lη2(η-1),N5=αη2(η-1)2+βη3(η-1)3

(1)

式(1)中N5為泡函數(shù)，α、β為任意常數(shù)，η=x/L，為單元的自然坐標(biāo)。

當(dāng)桿的內(nèi)力為壓力P時(shí)，由式(1)所示的形函數(shù)并利用虛功原理可得到平面剛架單元的變形剛度矩陣ke和幾何剛度矩陣ge，如式(2)的前兩式所示。

(2)

其中c1=EI/L3，c2=EA/L，a=4α2/5-12αβ/35+2β2/35，b=(α-3β/14)L，c=(44α2-22αβ+3β2)/77。令α=3，β=-3，則a=54/5，b=51L/14，c=621/77。

設(shè)平面剛架結(jié)構(gòu)的整體坐標(biāo)系為o-xy。x軸與x′軸間的夾角為φ。規(guī)定從x軸到x′ 軸以逆時(shí)針?lè)较驗(yàn)棣盏恼颉t局部坐標(biāo)系與整體坐標(biāo)系間的坐標(biāo)轉(zhuǎn)換矩陣T如式(2)的第三式所示。

對(duì)于空間剛架結(jié)構(gòu)，類似地采用一桿一單元，以桿件的一端為原點(diǎn)建立單元局部坐標(biāo)系o-x′y′z′，x′軸的正向指向桿的另一端，而y′軸和z′軸均為桿件橫截面的主慣性軸。桿件長(zhǎng)度為L(zhǎng)，橫截面積為A，彈性模量為E，剪切模量為G，橫截面對(duì)y′軸和z′軸的慣性矩為Iy和Iz，橫截面的極慣性矩為J。設(shè)空間剛架問(wèn)題的單元節(jié)點(diǎn)位移列向量為

q=(d1xd1yd1zφ1xφ1yφ1zd2xd2yd2zφ2xφ2yφ2z)T

(3)

根據(jù)文獻(xiàn)[3]中所示平面剛架剛度矩陣和空間剛架剛度矩陣的繼承與發(fā)展關(guān)系，并結(jié)合上述平面剛架泡函數(shù)有限元法的特點(diǎn)，得空間剛架單元的變形剛度矩陣ke和幾何剛度矩陣ge，如式(4)所示。

(4)

其中P為單元軸向壓力，a=EA/L，b=GJ/L，c=EIy/L3，d=EIz/L3，a1=4α2/5-12αβ/35+2β2/35，b1=(α-3β/14)L，c1=(44α2-22αβ+3β2)/77。令α=3，β=-3，則a1=54/5，b1=51L/14，c1=621/77。

設(shè)空間剛架結(jié)構(gòu)的整體坐標(biāo)系為o-xyz，nx′、ny′、nz′分別為各單元坐標(biāo)軸的方向余弦，則單元坐標(biāo)系與整體坐標(biāo)系間的坐標(biāo)轉(zhuǎn)換矩陣

(5)

根據(jù)以上的穩(wěn)定性理論并結(jié)合有限元法基本步驟，歸納出如下分岔點(diǎn)失穩(wěn)問(wèn)題的求解步驟。

(1)單位載荷作用下桿件內(nèi)力的計(jì)算

在研究點(diǎn)(將要施加待求載荷F)處施加與F同方向的單位載荷，根據(jù)有限元法的基本流程算得各桿件的內(nèi)力(P)，不必考慮結(jié)構(gòu)自重的影響。

(2)單元?jiǎng)偠染仃噆e和ge的計(jì)算與組裝

各單元都要計(jì)算ke，而只有當(dāng)P<0時(shí)才要計(jì)算ge(式(2)中的第2式以-P代替P)，接著將它們分別組裝成整體變形剛度矩陣K和整體幾何剛度矩陣G。

(3)考慮約束條件和求解廣義特征值問(wèn)題

剔除受約束的自由度，即從K和G中刪除該自由度所對(duì)應(yīng)的行和列，從而得到空間剛架結(jié)構(gòu)的變形方程為(K-FG)q=0，其具有非零解q的條件就是

|K-FG|=0

(6)

可利用逆迭代法來(lái)求解此廣義特征值問(wèn)題，得到最小特征值和相應(yīng)的特征向量。最小特征值即為臨界荷載Fcr，該特征向量即為屈曲模態(tài)，將此特征向量歸一化后再結(jié)合節(jié)點(diǎn)坐標(biāo)數(shù)據(jù)，可繪制屈曲模態(tài)實(shí)物圖。

(4)評(píng)估泡函數(shù)對(duì)計(jì)算結(jié)果精度的影響。

2 變幾何剛度法及失穩(wěn)判據(jù)

根據(jù)文獻(xiàn)[2]可知：和分岔點(diǎn)失穩(wěn)問(wèn)題處理法相比，在計(jì)算各單元的幾何剛度矩陣時(shí)，極值點(diǎn)失穩(wěn)問(wèn)題處理法也只考慮了待求載荷(F)，而忽略了其它載荷(如結(jié)構(gòu)自重、水平集中載荷、水平風(fēng)載荷等)的影響，不同點(diǎn)在于式(6)的右端為其它載荷的等效節(jié)點(diǎn)力向量(Fp)，被稱為結(jié)構(gòu)的屈曲測(cè)試向量。此時(shí)的空間剛架結(jié)構(gòu)的變形方程為

(K-FG)q=Fp

(7)

結(jié)構(gòu)的屈曲條件為q中至少有一個(gè)元素是無(wú)窮大量。為了增強(qiáng)可操作性，提出了如下失穩(wěn)判據(jù)：

|qi|max≥h/100,i=1…n

(8)

其中h表示剛架高度，n表示總體線位移自由度數(shù)。滿足此判據(jù)的最小F值即為Fcr。

因?yàn)槎家髊的系數(shù)矩陣行列式為零，所以上述兩種方法求得的臨界載荷值也相差不大。由于具有相同且固定的G矩陣，可將這兩種方法歸為一類，即定幾何剛度法。

結(jié)構(gòu)承受多種載荷作用時(shí)，宜采用極值點(diǎn)失穩(wěn)問(wèn)題處理法。當(dāng)待求載荷和其它載荷對(duì)結(jié)構(gòu)的變形具有相當(dāng)?shù)挠绊憰r(shí)，仍然采用處理極值點(diǎn)失穩(wěn)問(wèn)題的定幾何剛度法將會(huì)得出錯(cuò)誤的計(jì)算結(jié)果，給結(jié)構(gòu)設(shè)計(jì)帶來(lái)巨大的風(fēng)險(xiǎn)。為了克服這個(gè)缺點(diǎn)，處理極值點(diǎn)失穩(wěn)問(wèn)題的變幾何剛度法被提出。該法的變形方程為

(K-G)q=Fp

(9)

其中G矩陣包含了全部載荷的影響。變幾何剛度法和迭代法緊密結(jié)合，其具體實(shí)施步驟如下：

(1)采用分岔點(diǎn)失穩(wěn)問(wèn)題處理法計(jì)算出F的參考值F0，同時(shí)也生成K矩陣和Fp向量。

(2)采用變幾何剛度法計(jì)算出F的臨界值Fcr

取F0/2作為F的迭代初始值，上限值為3F0/2，步長(zhǎng)為單位力值?？紤]所有載荷的影響，計(jì)算出各桿件的內(nèi)力(P)，按式(4)計(jì)算單元幾何剛度矩陣，并組裝成總體幾何剛度矩陣G。式(9)中剔除受約束的自由度后，可解得總體節(jié)點(diǎn)位移向量q。這個(gè)迭代過(guò)程將循環(huán)地進(jìn)行下去，直至式(8)所示的判據(jù)得到滿足，此時(shí)所得的F值即為Fcr。

3 三點(diǎn)定位法

圖1為由角鋼搭建而成的空間剛架結(jié)構(gòu)[2]，確定各桿件單元坐標(biāo)系o-x′y′z′的方位是一個(gè)繁瑣而重要的問(wèn)題。圖2所示結(jié)構(gòu)為取自此類空間剛架的一部分。

圖1 鐵路控制塔計(jì)算模型圖

圖2 單元坐標(biāo)系示意圖

設(shè)在整體坐標(biāo)系中A點(diǎn)坐標(biāo)為(x1,y1,z1)，C點(diǎn)坐標(biāo)為(x2,y2,z2)，B點(diǎn)坐標(biāo)為(x3,y3,z3)。根據(jù)下面的一系列公式，可由A、C和B三點(diǎn)的整體坐標(biāo)確定斜桿AC的單元坐標(biāo)系。類似地，可用三點(diǎn)定位法確定斜桿DB和其它各斜桿的單元坐標(biāo)系。

(10)

(11)

4 算例

圖1中x-y平面為水平面，z方向?yàn)殂U垂方向。它具有28個(gè)連接點(diǎn)和96根桿件，其中24根立柱采用L50×5型角鋼(1#)，24根橫桿和48根斜桿均采用L56×5型角鋼(2#)。兩種角鋼截面的參數(shù)如表1所示。采用Q345型鋼材，其彈性模量為210 GPa，泊松比為0.3，密度為7 800 kg/m3。根據(jù)一桿一單元的網(wǎng)格劃分法可知，此剛架結(jié)構(gòu)具有28個(gè)節(jié)點(diǎn)和96個(gè)單元，25-28號(hào)節(jié)點(diǎn)固定。

(1)在1節(jié)點(diǎn)的-z軸方向上施加集中載荷F且不計(jì)結(jié)構(gòu)自重時(shí)，采用處理分岔點(diǎn)失穩(wěn)問(wèn)題的定幾何剛度法求F的臨界值Fcr1。

(2)在1節(jié)點(diǎn)的-z軸方向上施加集中載荷F且不計(jì)結(jié)構(gòu)自重時(shí)，采用處理極值點(diǎn)失穩(wěn)問(wèn)題的定幾何剛度法求F的臨界值Fcr2。

(3)在1節(jié)點(diǎn)的-z軸方向上施加集中載荷F和剛架背面(外法線方向?yàn)?x軸方向)承受水平風(fēng)壓載荷且考慮結(jié)構(gòu)自重時(shí)，采用處理極值點(diǎn)失穩(wěn)問(wèn)題的變幾何剛度法求F的臨界值Fcr3(單位長(zhǎng)度風(fēng)載荷p=-2ρbv2，風(fēng)速v=40 m/s，空氣密度ρ=1.29 kg/m3，b為角鋼斷面寬度)。

角鋼截面參數(shù)如表1所示，其中S表示角鋼截面號(hào)，q表示單位長(zhǎng)度的質(zhì)量，b為角鋼斷面寬度，A表示截面面積，J、Iy和Iz分別表示截面的扭轉(zhuǎn)慣性矩和兩個(gè)主軸慣性矩。節(jié)點(diǎn)坐標(biāo)數(shù)據(jù)如表2所示，其中N表示節(jié)點(diǎn)號(hào)，x,y和z表示節(jié)點(diǎn)的3個(gè)坐標(biāo)值。

表1 角鋼截面參數(shù)

表2 剛架的節(jié)點(diǎn)坐標(biāo)

對(duì)于24根立柱，有8種單元坐標(biāo)系；對(duì)于24根橫桿，有4種單元坐標(biāo)系。根據(jù)幾何關(guān)系可以直接寫出這12種單元坐標(biāo)系的方向余弦矩陣Λ。而對(duì)于48根斜桿，將三點(diǎn)定位算法編成規(guī)范的計(jì)算程序，可快速地獲得全部斜桿的Λ矩陣。

(1)分岔點(diǎn)失穩(wěn)(定幾何剛度法)

將上述泡函數(shù)有限元法和三點(diǎn)定位法理論編制成MATLAB程序，按前述分岔點(diǎn)失穩(wěn)問(wèn)題的求解步驟，最終求得此空間剛架結(jié)構(gòu)的臨界載荷Fcr1=951.87 kN。

當(dāng)采用普通有限元法即不考慮剛度矩陣(式(4))中的泡函數(shù)項(xiàng)時(shí)，適當(dāng)修改前述的MATLAB程序后求得此結(jié)構(gòu)的臨界載荷值為2 282.1 kN(一桿一單元)和952.57 kN(一桿二單元)。兩類算法計(jì)算結(jié)果的詳細(xì)比較如表3所示。由此可知：一桿一單元的泡函數(shù)有限元法可達(dá)到一桿二單元的普通有限元法的精度且計(jì)算時(shí)間更短；一桿一單元的普通有限元法的計(jì)算精度很差。

表3 兩類定幾何剛度分岔點(diǎn)失穩(wěn)有限元算法的對(duì)比

(2)極值點(diǎn)失穩(wěn)(定幾何剛度法)

根據(jù)式(7)可知K、G和Fp均不變，只需計(jì)算一次。采用迭代法求得臨界載荷Fcr2=951.88 kN。由此可知處理分岔點(diǎn)失穩(wěn)問(wèn)題的定幾何剛度法和處理極值點(diǎn)失穩(wěn)問(wèn)題的定幾何剛度法所求得的臨界載荷值相差不大，兩種方法沒(méi)有本質(zhì)區(qū)別。

(3)極值點(diǎn)失穩(wěn)(變幾何剛度法)

根據(jù)式(9)可知K和Fp均不變，只需計(jì)算一次，但G矩陣隨F的變化而變化。采用迭代法求得臨界載荷Fcr3=939.38 kN，明顯小于Fcr1和Fcr2。

由此可知，采用變幾何剛度法來(lái)處理極值點(diǎn)失穩(wěn)問(wèn)題具有安全性更高的優(yōu)點(diǎn)。

5 結(jié)論

空間剛架泡函數(shù)有限元法的突出優(yōu)點(diǎn)是既簡(jiǎn)潔又準(zhǔn)確。通過(guò)實(shí)際算例，如下的三大創(chuàng)新點(diǎn)得到了驗(yàn)證。

(1)采用泡函數(shù)法，在保證計(jì)算精度的前提下大幅降低了問(wèn)題的自由度數(shù)；

(2)提出變幾何剛度法和失穩(wěn)判據(jù)，全面考慮了各種載荷因素對(duì)幾何剛度矩陣的影響，優(yōu)化了極值點(diǎn)失穩(wěn)理論；

(3)采用三點(diǎn)定位法，簡(jiǎn)化了單元坐標(biāo)系的定位問(wèn)題。

而目前的各種商用有限元分析軟件都不具有上述功能，因此本文所論述的空間剛架穩(wěn)定性分析的泡函數(shù)有限元法具有極大的實(shí)用價(jià)值。

[1] 高軒能,王仁輝.框架結(jié)構(gòu)穩(wěn)定分析的泡函數(shù)有限元法[J].華僑大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,27(2):155-158.

[2] 龍馭球,包世華,匡文起,等.結(jié)構(gòu)力學(xué)Ⅱ(第2版)[M].北京:高等教育出版社,2006.

[3] DARYL L.LOGAN.有限元方法基礎(chǔ)教程(第五版)[M].張榮華,王藍(lán)婧,李繼榮,等譯,北京:電子工業(yè)出版社,2014.

[4] 劉文國(guó),任文敏,張維.用復(fù)合有限條-氣泡函數(shù)法研究加肋殼體穩(wěn)定性[J].清華大學(xué)學(xué)報(bào)(自然科學(xué)版),1999,39(11):61-64.

[5] 夏擁軍,陸念力.梁桿結(jié)構(gòu)穩(wěn)定性分析的高精度Euler-Bernoulli梁?jiǎn)卧猍J].沈陽(yáng)建筑大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,22(3):362-366.

[6] 劉慧娟,羅永峰,楊綠峰,等.單層網(wǎng)殼結(jié)構(gòu)穩(wěn)定性分析的隨機(jī)缺陷模態(tài)疊加法[J].同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版),2012,40(9):1294-1299.

[7] 周 飛,宛樹(shù)旗.單層網(wǎng)殼結(jié)構(gòu)穩(wěn)定性分析[J].國(guó)外建材科技,2007,28(3):119-122.

[8] 彭細(xì)榮,楊慶生,孫卓.有限單元法及其應(yīng)用[M].北京:清華大學(xué)出版社,2012.

[9] 徐斌,高躍飛,余龍.MATLAB有限元結(jié)構(gòu)動(dòng)力學(xué)分析與工程應(yīng)用[M].北京:清華大學(xué)出版社,2009.

[10]黃志新,劉成柱.ANSYS Workbench 14.0超級(jí)學(xué)習(xí)手冊(cè)[M].北京:人民郵電出版社,2014.

[11]龍述堯,蒯形成,劉騰喜.計(jì)算力學(xué)[M].長(zhǎng)沙:湖南大學(xué)出版社,2007.

[12]P I KATTAN.MATLAB有限元分析與應(yīng)用[M].韓來(lái)彬譯,北京:清華大學(xué)出版社,2006.

[13]CHAJES A.Principles of Structural Stability Theory[M].New York:Prentice-Hall Inc,1974.

[14]WEAVER W,JOHNSTON P R.Finite element for ftructural analysis[M].New York:Prentice-Hall Inc,1984.

[15]CAO JIWEI,LIU MINGFANG,CHEN SHAOCHUN.Mixed Finite Element Formats of any Order Based on Bubble Functions for Stationary Stokes Problem[J].Journal of Math,2016,31(1):87-95.

Abubble-functionfiniteelementmethodforthestabilityanalysisofspaceframestructures

CHEN Liang，CHEN Zhen-zhong
(Faculty of Aerospace Engineering,Shenyang Aerospace University,Shenyang 110136,China)

Stability of space frame structures was analyzed by a Bubble-Function Finite Element Method(BFFEM)in this paper.Both a variable geometric stiffness method and a buckling criterion were presented with the consideration of various load influences on the geometric stiffness matrix.A formal program was written by a Three Points Method to define different element coordinate systems.For one bar one element,a very precise critical load for the space frame was calculated using the above three theories,which indicated that the BFFEM was simple and accurate for analyzing the stability of space frame structures.

space frame structures;stability analysis;bubble-function;variable geometric stiffness;three points method;finite element method

2017-10-08

陳 梁(1988-)，男，湖南長(zhǎng)沙人，碩士研究生，主要研究方向：結(jié)構(gòu)動(dòng)力學(xué)分析與有限元法應(yīng)用，E-mail:chenliang0102@163.com；陳振中(1963-)，男，遼寧沈陽(yáng)人，教授，博士，主要研究方向：結(jié)構(gòu)疲勞斷裂與可靠性分析，E-mail:zhenzhong_chen@hotmail.com。

2095-1248(2017)06-0040-06

O302;TU328

A

10.3969/j.issn.2095-1248.2017.06.007

吳萍 英文審校：趙歡)