基于CR列式的動(dòng)力非線性有限元程序開發(fā)

王 濤，沈銳利

(西南交通大學(xué) 土木工程學(xué)院，四川 成都 610031)

?

基于CR列式的動(dòng)力非線性有限元程序開發(fā)

王 濤，沈銳利

(西南交通大學(xué) 土木工程學(xué)院，四川 成都 610031)

基于CR列式非線性有限元計(jì)算理論，提出了針對桿、梁結(jié)構(gòu)的Newmark-β非線性有限元?jiǎng)恿r(shí)程積分算法,在動(dòng)力時(shí)程計(jì)算中，可以考慮結(jié)構(gòu)的幾何非線性效應(yīng)，通過計(jì)算單元伸長并扣除結(jié)構(gòu)剛體位移得到單元的實(shí)際變形和內(nèi)力；闡述了非線性有限元?jiǎng)恿r(shí)程計(jì)算程序的計(jì)算原理,開發(fā)了計(jì)算程序，建立有限元模型；通過算例與ANSYS的計(jì)算結(jié)果進(jìn)行了對比。結(jié)果表明：程序非線性靜力與動(dòng)力時(shí)程計(jì)算結(jié)果與ANSYS差別很小，且計(jì)算速度更快，程序算法正確可靠，適用于桿、梁結(jié)構(gòu)大幅度非線性振動(dòng)的計(jì)算。

橋梁工程；有限元程序；CR列式；幾何非線性；動(dòng)力時(shí)程；斜拉索

0 引 言

考慮結(jié)構(gòu)幾何非線性有限元計(jì)算方法是工程分析中的熱點(diǎn)。對幾何非線性，國內(nèi)外學(xué)者提出了多種方法。最常用的是以結(jié)構(gòu)變形前為參考，建立平衡方程的全拉格朗日法(TL列式法)和以結(jié)構(gòu)變形后為參考，建立平衡方程的更新拉格朗日法(UL列式法)。K.J.Bathe等[1-2]在大量研究基礎(chǔ)上，建立了三維梁單元大位移、大轉(zhuǎn)動(dòng)、小應(yīng)變的UL列式和TL列式分析方法。陳政清等[3]改進(jìn)了K.J.Bathe建立的非線性梁單元，提高了UL列式計(jì)算效率。

對于大轉(zhuǎn)動(dòng)問題，為了得到更為精確的結(jié)果，使用TL或UL列式往往需要增加荷載步，從而顯著增加了計(jì)算時(shí)間，且因涉及到大轉(zhuǎn)動(dòng)問題，UL列式方法的推導(dǎo)過程往往都非常復(fù)雜，在實(shí)際編程的應(yīng)用中難度較大。國內(nèi)外廣大學(xué)者對該問題其它解決途徑進(jìn)行了研究，流動(dòng)坐標(biāo)系方法——CR列式法(Co-rotational Formulation)成為了結(jié)構(gòu)幾何非線性有限元計(jì)算另一關(guān)注點(diǎn)。

G.A.Wempner等[4]，T.Belytschko等[5]分別提出了CR列式的概念。A.Crisfield等[6]對不同單元幾何非線性CR列式提出了一致列式方法。K.Hsiao等[7]，C.A.Felippa等[8]對梁單元的CR列式計(jì)算方法進(jìn)行了深入研究。周凌遠(yuǎn)等[9]將UL列式計(jì)算概念與CR列式方法結(jié)合起來，提出了基于UL法的CR列式三維梁單元計(jì)算方法。潘永仁[10]，唐茂林[11]基于CR列式幾何非線性計(jì)算方法開發(fā)了有限元程序，運(yùn)用到大跨度懸索橋施工過程監(jiān)控，取得了良好的結(jié)果。

以上研究主要關(guān)注結(jié)構(gòu)靜力計(jì)算幾何非線性問題。目前，在結(jié)構(gòu)動(dòng)力時(shí)程計(jì)算中，通常使用振型分解法或線性Newmark-β法來得到結(jié)構(gòu)的振動(dòng)響應(yīng)，對于結(jié)構(gòu)的小幅度振動(dòng)，線性分析方法尚能達(dá)到滿足工程要求的精度，但當(dāng)大跨度或柔性結(jié)構(gòu)發(fā)生較大幅度振動(dòng)時(shí)(特別是索結(jié)構(gòu))，振動(dòng)往往是非線性的，使用線性的計(jì)算方法得到的結(jié)果可能與實(shí)際情況偏差較大。

為了在動(dòng)力時(shí)程計(jì)算中模擬結(jié)構(gòu)的動(dòng)態(tài)受力情況，筆者基于CR列式對單元內(nèi)力的計(jì)算方法，建立了考慮結(jié)構(gòu)幾何非線性的Newmark-β動(dòng)力時(shí)程計(jì)算方法，編制了計(jì)算程序，使用ANSYS驗(yàn)證了程序的正確性與可靠性。相對于一般基于UL列式幾何非線性計(jì)算方法的商業(yè)通用有限元軟件(如ANSYS等)，該方法計(jì)算速度快，力學(xué)概念簡潔(通常不必計(jì)算、推導(dǎo)復(fù)雜的單元非線性大位移剛度矩陣)，適合于模擬柔性結(jié)構(gòu)(特別是索結(jié)構(gòu))的大幅度非線性振動(dòng)。

1 CR列式方法基本原理

以2維梁單元為例，在圖1(a)中給出了在結(jié)構(gòu)總體坐標(biāo)系XY中未變形的二維梁單元，單元受到總體坐標(biāo)系下的節(jié)點(diǎn)力以后，桿端的位移與轉(zhuǎn)角如圖1(b)。

圖1 基于CR列式方法的單元計(jì)算示意

根據(jù)CR列式計(jì)算原理，在單元上附加一個(gè)流動(dòng)坐標(biāo)系X1Y1，流動(dòng)坐標(biāo)系原點(diǎn)始終位于單元i節(jié)點(diǎn)上，X1軸始終沿單元節(jié)點(diǎn)ij方向，Y1軸始終垂直于X1軸。

單元變形后，流動(dòng)坐標(biāo)系X1Y1位置與單元在流動(dòng)坐標(biāo)系內(nèi)的伸長與轉(zhuǎn)角如圖1(b)。

對于圖1(a)可以計(jì)算得到單元的長度與傾角：

(1)

對于圖1(b)可以計(jì)算得到單元的長度與傾角：

(2)

對于梁單元，在單元流動(dòng)坐標(biāo)X1Y1中，單元節(jié)點(diǎn)的實(shí)際位移為單元伸長后的長度與初始長度的差。單元節(jié)點(diǎn)的實(shí)際轉(zhuǎn)角為總體坐標(biāo)系下節(jié)點(diǎn)的轉(zhuǎn)角值減去單元在整體坐標(biāo)系下發(fā)生的剛體轉(zhuǎn)動(dòng)值，顯然，圖1中剛體轉(zhuǎn)動(dòng)值為α2-α1。

對于平面梁單元，筆者使用了較為簡單的剛體轉(zhuǎn)動(dòng)計(jì)算公式。剛體轉(zhuǎn)動(dòng)也可以根據(jù)單元變形前后X1軸基矢量的指向，使用矢量法則來計(jì)算[9-10]。

綜上所述，如圖1(b)，在流動(dòng)坐標(biāo)系中，單元變形后的真實(shí)伸長與轉(zhuǎn)角具體表達(dá)式為：

(3)

每一步迭代中梁單元在單元隨動(dòng)坐標(biāo)系X1Y1中的節(jié)點(diǎn)的實(shí)際位移向量為：

{u}e={0 0 θi1uj10 θj1}

(4)

式(4)中，不計(jì)轉(zhuǎn)動(dòng)項(xiàng)，即為直桿單元的位移向量。由式(3)、式(4)可以看出，對于直桿單元，由于在流動(dòng)坐標(biāo)系中僅有X1軸方向的變形，使用CR列式方法計(jì)算得到的單元伸長在幾何上是精確的。對于梁單元，得到的單元伸長，通過扣除剛體轉(zhuǎn)動(dòng)得到節(jié)點(diǎn)在單元坐標(biāo)系內(nèi)的轉(zhuǎn)動(dòng)值也是幾何精確的?？梢缘玫骄植孔鴺?biāo)系下單元實(shí)際內(nèi)力：

{f}e=[k]e{u}e

(5)

式中：{f}e為單元內(nèi)力的節(jié)點(diǎn)力向量；[k]e為單元彈性剛度矩陣。

一般情況下考慮結(jié)構(gòu)幾何非線性計(jì)算時(shí)，可認(rèn)為整體結(jié)構(gòu)處于大變形、大轉(zhuǎn)動(dòng)、小應(yīng)變狀態(tài)，那么在單元局部流動(dòng)坐標(biāo)系內(nèi)可認(rèn)為力與位移的關(guān)系是線性的，[k]e可取為線彈性的桿、梁單元矩陣[10-11]。

在非線性計(jì)算中，總體結(jié)構(gòu)在外力作用下變形后的總體節(jié)點(diǎn)內(nèi)力向量為位移的非線性函數(shù)，可表示為{R({u})}，總體節(jié)點(diǎn)內(nèi)力向量流程如圖2。

圖2 有限元模型總體節(jié)點(diǎn)內(nèi)力向量計(jì)算流程

圖3 有限元模型總體節(jié)點(diǎn)外力向量計(jì)算流程

根據(jù)圖2與圖3的計(jì)算流程，可以認(rèn)為，在一般的CR列式算法中，結(jié)構(gòu)在總體坐標(biāo)系下，所有的幾何非線性效應(yīng)都是由位移造成的坐標(biāo)轉(zhuǎn)換矩陣的變化來反映的。由于CR列式為幾何精確方法，所以從理論上來講，荷載的分級(jí)加載的步數(shù)對最后結(jié)果無影響，差別來自于數(shù)值舍入的誤差。

如果結(jié)構(gòu)處于最終的靜力平衡狀態(tài)，則有：

(6)

對于式(6)可以使用迭代計(jì)算得到結(jié)構(gòu)在外力作用下達(dá)到靜力平衡狀態(tài)時(shí)總體坐標(biāo)系下各個(gè)節(jié)點(diǎn)的位移向量{u}。關(guān)于桿系CR列式非線性有限元靜力計(jì)算方法的細(xì)節(jié)可以參考文獻(xiàn)[9-11]。

2 基于CR列式的非線性有限元?jiǎng)恿r(shí)程積分

2.1 基本理論

筆者開發(fā)了基于CR列式的非線性Newmark-β有限元?jiǎng)恿r(shí)程積分方法，其基本假定與普通的Newmark-β法相同，可以寫為：

(7)

(8)

其中各個(gè)參數(shù)具體表達(dá)式可以參見文獻(xiàn)[12]。考慮幾何非線性時(shí)，結(jié)構(gòu)t+Δt時(shí)刻的振動(dòng)方程為：

(9)

將式(7)，式(8)代入式(9)，可以得到：

(10)

這樣，有限元非線性動(dòng)力時(shí)程積分就轉(zhuǎn)化為在每一個(gè)時(shí)間步內(nèi)求解非線性方程組(10)的問題。

式(10)左端中的總體節(jié)點(diǎn)內(nèi)力向量{R}是總體節(jié)點(diǎn)位移向量{u}的非線性函數(shù)。對于{R}可以按照圖 2中的流程來計(jì)算 。

在考慮結(jié)構(gòu)幾何非線性時(shí)，式(10)右端中，{F}表示結(jié)構(gòu)承受的總體節(jié)點(diǎn)外力向量，{Fe}表示結(jié)構(gòu)初始單元力導(dǎo)致的等效總體節(jié)點(diǎn)外力向量(如ANSYS中初應(yīng)變的概念)，它們都是位移{u}的非線性函數(shù)，可以按照圖3中的流程來計(jì)算。

式(10)左右兩端還存在與總體質(zhì)量矩陣[M] 、總體阻尼矩陣[C] 相關(guān)的等效力項(xiàng)。在非線性動(dòng)力時(shí)程計(jì)算中，節(jié)點(diǎn)位置的變化必然會(huì)導(dǎo)致[M]的變化，所以在每一個(gè)時(shí)間步的計(jì)算中要根據(jù)當(dāng)前節(jié)點(diǎn)位置使用坐標(biāo)轉(zhuǎn)換矩陣來更新[M]。對于阻尼，由于其復(fù)雜性，計(jì)算中不考慮[C]的變化。

可以使用牛頓迭代法來求解式(10)，當(dāng)?shù)词諗繒r(shí)，式(10)左右兩端是不相等的，計(jì)算兩端的差值就可以得到迭代的不平衡力。在迭代過程中使用的切線剛度矩陣為Newmark-β法等效總體剛度矩陣：

[K]=a0[M]+a1[C]+[Kk]

(11)

式中：[M],[C],[Kk]分別為將節(jié)點(diǎn)坐標(biāo)更新到當(dāng)前迭代步位置上時(shí)，結(jié)構(gòu)的總體質(zhì)量矩陣、總體阻尼矩陣、線性靜力計(jì)算時(shí)計(jì)算的總體剛度矩陣。在計(jì)算中[Kk] 使用單元彈性剛度矩陣疊加單元應(yīng)力剛度矩陣來組集。為了加快迭代計(jì)算收斂速度，式(11)也根據(jù)節(jié)點(diǎn)位置狀態(tài)來更新。

2.2 單元計(jì)算

筆者使用的梁單元彈性剛度矩陣為：

(12)

式中：l0為單元的無應(yīng)力長度；E為彈性模量；A為單元截面積；Iz為單元抗彎慣性矩。

若Iz=0則式(12)為平面桿單元的彈性剛度矩陣。若單元上存在軸力，則會(huì)產(chǎn)生應(yīng)力剛度矩陣，在CR列式非線性計(jì)算迭代中，應(yīng)力剛度矩陣一般起到加快收斂速度的作用，對最后的計(jì)算結(jié)果無影響[11]。梁單元應(yīng)力剛度矩陣為：

(13)

式中：f為單元軸力。

(14)

式中：l為單元當(dāng)前狀態(tài)下的長度；ρ 為單元材料的質(zhì)量密度。

桿單元應(yīng)力剛度矩陣為：

(15)

桿單元一致質(zhì)量矩陣為：

(16)

桿、梁單元的坐標(biāo)轉(zhuǎn)換矩陣為：

(17)

CR列式的非線性計(jì)算中，單元坐標(biāo)轉(zhuǎn)換矩陣是隨著流動(dòng)坐標(biāo)系的變化而更新的。根據(jù)每一步迭代得到的位移值更新有限元模型節(jié)點(diǎn)坐標(biāo)，再根據(jù)每個(gè)單元節(jié)點(diǎn)i, j的位置與單元長度l可計(jì)算得到：

(18)

2.3 算法流程

綜上所述，筆者基于CR列式的幾何非線性Newmark-β法有限元?jiǎng)恿r(shí)程積分計(jì)算流程如圖4。

圖4 非線性有限元?jiǎng)恿r(shí)程計(jì)算流程

為了更清晰地闡述筆者編程原理，使用了較為簡單的平面桿、梁單元來構(gòu)建程序。如需要使程序具有更廣泛的適用性，可將程序擴(kuò)展為三維空間非線性動(dòng)力時(shí)程計(jì)算程序，編程思路與圖4相同。

如果單元為三維空間桿單元，根據(jù)CR列式原理，只需要計(jì)算單元在三維空間的伸長便可以獲得單元的實(shí)際內(nèi)力。但對于梁、殼單元，三維空間中剛體轉(zhuǎn)動(dòng)情況是較為復(fù)雜的，為了計(jì)算三維單元扣除剛體轉(zhuǎn)動(dòng)的真實(shí)節(jié)點(diǎn)位移獲得單元實(shí)際內(nèi)力，通常要使用三維矢量代數(shù)公式來計(jì)算，文獻(xiàn)[6-10]中較完整地討論了這個(gè)問題 的技術(shù)解決方案。

筆者只闡述了使用的關(guān)鍵技術(shù)方法，更多的關(guān)于有限元程序開發(fā)技術(shù)，如：總體剛度矩陣組裝、節(jié)點(diǎn)約束、節(jié)點(diǎn)自由度釋放、節(jié)點(diǎn)耦合、矩陣求解、收斂條件判斷、節(jié)點(diǎn)強(qiáng)制位移等可以參考文獻(xiàn)[12]或其他關(guān)于有限元方法的著作。

3 算例驗(yàn)證

筆者使用MATLAB數(shù)值計(jì)算平臺(tái)，開發(fā)了基于CR列式的桿系結(jié)構(gòu)非線性動(dòng)力有限元程序，程序特點(diǎn)如下：①可以考慮結(jié)構(gòu)幾何非線性計(jì)算結(jié)構(gòu)的靜力狀態(tài)；②可以考慮結(jié)構(gòu)大變形后的靜力狀態(tài)，計(jì)算結(jié)構(gòu)的動(dòng)力特性；③可以計(jì)算結(jié)構(gòu)的非線性動(dòng)力時(shí)程；④可以根據(jù)非線性動(dòng)力時(shí)程計(jì)算結(jié)果繪制輸出結(jié)構(gòu)非線性振動(dòng)的動(dòng)畫圖形。

為了驗(yàn)證筆者程序的正確性與可靠性，將算例的結(jié)果與ANSYS計(jì)算結(jié)果進(jìn)行對比。程序計(jì)算中的重力加速度均取為g=9.8 m/s2。計(jì)算中不設(shè)置結(jié)構(gòu)阻尼，由文獻(xiàn)[13]可知ANSYS計(jì)算中默認(rèn)的算法阻尼不為0，即：Newmark-β法中γ=0.505，所以在本文程序計(jì)算中使用相同設(shè)置。ANSYS中分別使用Link1與 Beam3單元模擬平面桿、梁結(jié)構(gòu)。

3.1 算例1

工程中常見的索-梁組合結(jié)構(gòu)有限元模型如圖5，主梁分為4個(gè)梁單元，拉索共分為10個(gè)直桿單元，彈性模量E、材料質(zhì)量密度ρ，單元截面積A，梁單元抗彎剛度Iz，桿單元初始軸力H，單元物理參數(shù)如表1。

圖5 索-梁組合結(jié)構(gòu)有限元模型

表1 單元物理參數(shù)

考慮幾何非線性計(jì)算結(jié)構(gòu)在自重作用下的靜力構(gòu)型，得到ANSYS與筆者計(jì)算結(jié)果對比，如圖6。

圖6 索-梁組合結(jié)構(gòu)靜力計(jì)算ANSYS與本文程序計(jì)算結(jié)果

從圖6可以看出，本文非線性靜力計(jì)算結(jié)果與ANSYS差別很小。計(jì)算靜力構(gòu)型下的動(dòng)力特性，得到結(jié)構(gòu)前2階模態(tài)如圖7。

圖7 結(jié)構(gòu)前2階模態(tài)

靜力計(jì)算結(jié)果得到拉索的平均軸力為50.53 kN，由文獻(xiàn)[14]的計(jì)算公式可以得到拉索局部的1階自振頻率為3.173 Hz。由圖7及文獻(xiàn)[14]可知，由于拉索的1階自振頻率接近整體結(jié)構(gòu)的1階自振頻率，整體結(jié)構(gòu)發(fā)生振動(dòng)時(shí)，拉索會(huì)在端點(diǎn)位移激勵(lì)(節(jié)點(diǎn)5)在拉索局部坐標(biāo)系Y1方向的分量作用下發(fā)生1階非線性主共振，文獻(xiàn)[14]稱之為索-梁相關(guān)振動(dòng)。

在節(jié)點(diǎn)5上作用Y方向豎向力P=-10 kN，靜力計(jì)算后釋放節(jié)點(diǎn)力，動(dòng)力時(shí)程積分取時(shí)間步長為0.02 s，分別使用ANSYS與本文程序計(jì)算結(jié)構(gòu)在有自重狀態(tài)下的振動(dòng)時(shí)程[13]，如圖8。

圖8 索-梁組合結(jié)構(gòu)非線性振動(dòng)ANSYS與本文程序計(jì)算結(jié)果

從圖8可以看出，使用ANSYS計(jì)算的結(jié)果與本文計(jì)算得到的振動(dòng)時(shí)程圖是接近的。在釋放節(jié)點(diǎn)力后結(jié)構(gòu)發(fā)生了振動(dòng)，索-梁相關(guān)振動(dòng)[14]導(dǎo)致了拉索的共振，拉索端部約0.01 m幅值的位移激勵(lì)，導(dǎo)致拉索中點(diǎn)發(fā)生了0.06 m的振動(dòng)。結(jié)構(gòu)振動(dòng)體現(xiàn)了“拍振”的非線性振動(dòng)性質(zhì)[15]，振動(dòng)能量在拉索與整體結(jié)構(gòu)之間互相傳遞，呈“此消彼長”的趨勢。

如果程序存在錯(cuò)誤，那么隨著時(shí)間的發(fā)展，振動(dòng)時(shí)程圖的差別會(huì)變得更為明顯。所以提取拉索1/2點(diǎn)(節(jié)點(diǎn)11)在30～40 s之間Y1方向振動(dòng)時(shí)程對比ANSYS與本文程序計(jì)算結(jié)果，如圖9。

圖9 ANSYS與本文程序結(jié)果對比

從圖9可以看出，對于30～40s時(shí)間段內(nèi)節(jié)點(diǎn)11的振動(dòng)時(shí)程，由于算法上的差別，筆者計(jì)算得到的振動(dòng)幅值較ANSYS略微偏大，但兩者的結(jié)果總體是很接近的，這說明了筆者程序的正確性。

3.2 算例2

某實(shí)際斜拉橋拉索(圖10)的兩端點(diǎn)直線距離l=112.775 m，傾角θ=24°，每延米質(zhì)量m=60.1 kg/m，彈性模量E=1.96×1011Pa， 拉索初始軸力H=3 665 kN，截面積A=0.007 5 m2，拉索共分為20個(gè)直桿單元。

圖10 斜拉索在端點(diǎn)位移激勵(lì)下的有限元模型

考慮幾何非線性計(jì)算拉索在自重作用下的靜力構(gòu)型，得到ANSYS與本文計(jì)算結(jié)果對比，如圖11。

圖11 拉索靜力計(jì)算ANSYS與本文程序計(jì)算結(jié)果

從圖11可以看出，本文非線性靜力計(jì)算結(jié)果與ANSYS差別很小。計(jì)算拉索在靜力構(gòu)型下的動(dòng)力特性，得到結(jié)構(gòu)前2階模態(tài)，如圖12。

圖12 拉索前2階模態(tài)

如圖10，設(shè)拉索端點(diǎn)位移激勵(lì)U=ΔUsin(ωt) ，取端點(diǎn)位移激勵(lì)幅值為ΔU=0.1 m，使用強(qiáng)制位移施加端點(diǎn)位移激勵(lì)，激勵(lì)頻率ω為1.105×2=2.210 Hz ，動(dòng)力時(shí)程積分時(shí)間步長取為0.02 s。

由于端點(diǎn)位移激勵(lì)接近拉索的1階自振頻率的2倍，拉索會(huì)在端點(diǎn)位移激勵(lì)X1方向分量作用下發(fā)生1階參數(shù)主共振，拉索垂度較小，2階頻率接近1階頻率的2倍，拉索會(huì)在端點(diǎn)激勵(lì)Y1方向分量作用下發(fā)生2階主共振[15]。得到分別使用ANSYS與筆者程序計(jì)算拉索在有重力狀態(tài)下1/4點(diǎn)(節(jié)點(diǎn)6)、1/2點(diǎn)(節(jié)點(diǎn)11)沿Y1方向振動(dòng)時(shí)程，如圖13。

圖13 拉索非線性振動(dòng)ANSYS與本文程序計(jì)算結(jié)果

由圖13可以看出，由于拉索上發(fā)生了2階主共振與1階參數(shù)共振，拉索1/4點(diǎn)與1/2點(diǎn)都發(fā)生了較大幅度的振動(dòng)。提取90～100 s內(nèi)拉索1/2點(diǎn)Y1方向振動(dòng)時(shí)程，對比ANSYS與本文程序計(jì)算結(jié)果，如圖14。

圖14 ANSYS與本文程序結(jié)果對比

由圖14可知，本文程序計(jì)算結(jié)果與ANSYS差別很小，驗(yàn)證了程序的正確性。

筆者程序可以根據(jù)計(jì)算結(jié)果生成振動(dòng)時(shí)程的動(dòng)畫，可以更為直觀地觀察拉索非線性振動(dòng)的發(fā)展過程，但限于表達(dá)方式，圖15列出了拉索在各個(gè)時(shí)刻非線性振動(dòng)的振動(dòng)形狀。

圖15 拉索各個(gè)時(shí)間點(diǎn)振動(dòng)形狀(有放大)

由圖15可以看出，拉索在端點(diǎn)位移激勵(lì)作用下，首先發(fā)生了2階主共振，隨著時(shí)間的發(fā)展，拉索1/2點(diǎn)的振幅變大，拉索發(fā)生了1階參數(shù)共振，這兩種大幅度非線性振動(dòng)在拉索上是同時(shí)存在的。

由于算例1、算例2均是計(jì)算結(jié)構(gòu)在自重靜力構(gòu)型下的振動(dòng)響應(yīng)，所以結(jié)構(gòu)靜力變形后的狀態(tài)為振動(dòng)平衡位置，振動(dòng)時(shí)程的初始位移不為0。

4 結(jié) 論

1)筆者提出了基于CR列式的非線性動(dòng)力時(shí)程有限元算法，詳細(xì)地闡述了程序開發(fā)過程。編制了有限元計(jì)算程序，通過算例的靜力與動(dòng)力計(jì)算與ANSYS結(jié)果的對比，驗(yàn)證了程序的正確、可靠性。

2)相對于ANSYS這樣使用UL列式計(jì)算結(jié)構(gòu)內(nèi)力的通用軟件，基于筆者算法開發(fā)的有限元程序具有更加簡潔、實(shí)用、高效的特點(diǎn)。在算例1動(dòng)力時(shí)程計(jì)算中，計(jì)算2 000步，使用ANSYS耗時(shí)約63 s，筆者程序耗時(shí)約18 s。在算例2動(dòng)力時(shí)程計(jì)算中，由于使用了強(qiáng)制位移加載，計(jì)算5 000步ANSYS耗時(shí)約為430 s，筆者程序耗時(shí)約50 s。

3)本文程序可以計(jì)算在外力作用下柔性桿系結(jié)構(gòu)的非線性振動(dòng)，也可以計(jì)算結(jié)構(gòu)在強(qiáng)制位移激勵(lì)作用下的非線性振動(dòng)。非線性有限元?jiǎng)恿r(shí)程計(jì)算更好地反映了結(jié)構(gòu)的受力細(xì)節(jié)。筆者方法為桿系結(jié)構(gòu)非線性振動(dòng)研究提供了較好的編程解決方案，適用于柔性結(jié)構(gòu)，特別是索結(jié)構(gòu)非線性振動(dòng)的研究。

[1] Bathe K J,Ramm E,Wilson E L.Finite element formulations for large deformation dynamic analysis [J].International Journal for Numerical Methods and Engineering,1975,9(2):353-386.

[2] Bathe K J,Bolourchi S.Large displacement analysis of three-dimensional beam structures [J].International Journal for Numerical Methods and Engineering,1979,14(7):961-986.

[3] 陳政清,曾慶元,顏全勝.空間桿系結(jié)構(gòu)大撓度問題內(nèi)力分析的UL列式法[J].土木工程學(xué)報(bào),1992,25(5):34-44. Chen Zhengqing,Zeng Qingyuan,Yan Quansheng.A UL formulation for internal force analysis of spatial frame structures with large displacement [J].China Civil Engineering Journal,1992,25(5):34-44.

[4] Wempner G A.Finite elements,finite rotations and small strains of flexible shells [J].International Journal of Solids and Structures,1969,5(7):117-153.

[5] Belytschko T,Hsieh B J.Non-linear transient finite element analysis with convected co-ordinates [J].International Journal for Numerical Methods and Engineering,1973,7(9):255-271.

[6] Crisfieldm A,Moita G F.A unified co-rotational framework for solids,shells and beams [J].International Journal of Solids and Structures,1996,33(24):2969-2992.

[7] Hsiao K,Linw Y.Co-rotational finite element formulation for buckling and post buckling analyses of spatial beams [J].Computer Methods in Applied Mechanics and Engineering,2000,188(3):567-594.

[8] Felippa C A,Haugen B.A unified formulation of small-strain co-rotational finite elements[ J].Theory Computer Methods in Applied Mechanics and Engineering,2005,194(19):2285-2 335.

[9] 周凌遠(yuǎn),李喬.基于UL法的CR列式三維梁單元計(jì)算方法[J].西南交通大學(xué)學(xué)報(bào),2006,41(6):690-695. Zhou Lingyuan,Li Qiao.Updated lagrangian co-rotational formulation for geometrically nonlinear FE analysis of 3-D beam element [J].Journal of Southwest Jiaotong University,2006,41(6):690-695.

[10] 潘永仁.懸索橋結(jié)構(gòu)非線性分析理論與方法[M].北京：人民交通出版社,2001. Pan Yongren.Non-Linear Analysis Theory Method for Suspension Bridge Structure [M].Beijing:China Communications Press,2001.

[11] 唐茂林.大跨度懸索橋空間幾何非線性分析與軟件開發(fā)[D].成都：西南交通大學(xué),2003. Tang Maolin.3D Geometric Nonlinear Analysis of Long-Span Suspension bridge and its Software Development [D].Chengdu:Southwest Jiaotong University,2003.

[12] 徐榮橋.結(jié)構(gòu)分析的有限元法與MATLAB程序設(shè)計(jì)[M].北京：人民交通出版社,2006. Xu Rongqiao.Finite Element Method in Structural Analysis and MATLAB Programming [M].Beijing:China Communications Press,2006.

[13] 王新敏.ANSYS結(jié)構(gòu)動(dòng)力分析與應(yīng)用[M].北京：人民交通出版社.2014. Wang Xinmin.Structural Dynamic Analysis and Application with ANSYS [M].Beijing:China Communications Press.2014.

[14] 王濤,沈銳利.斜拉橋索-梁相關(guān)振動(dòng)概念與研究方法初探[J].振動(dòng)與沖擊,2013,20(32):29-34. Wang Tao,Shen Ruili.Primary investigation on the concept and method of cable-beam vibration in cable-stayed bridge [J].Journal of Vibration and Shock,2013,20(32):29-34.

[15] Nayfeh A.H,Mook D T.Nonlinear Oscillations [M].New York:Wiley Press,1984.

Developing of Dynamic Nonlinear Finite Element Method ProgramBased on Co-Rotational Formulation

Wang Tao, Shen Ruili

(School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, Sichuan, China)

Based on the nonlinear FEM computational theory of co-rotational (CR) formulation, a Newmark-βnonlinear FEM dynamic time-history algorithm for the structures of truss and beam was proposed. The geometric nonlinearity of the structures could be considered in the time-history calculation. The real displacement and internal force of the elements could be obtained by calculating the element extension and by deducting the rigid body displacement of the structure. The calculation principles of the geometric nonlinearity dynamic program were elaborated. A calculation program was developed and the FEM model was established. The comparison was carried out between the example and the calculation results of the ANSYS. The results show that the difference of the calculation results between the program and ANSYS on nonlinear static and dynamic time-history is very small, and the program is faster than ANSYS, which is applicative for large amplitude vibration calculation of truss and beam structures.

bridge engineering; FEM program; co-rotational formulation; geometric nonlinearity; dynamic time-history; stayed cable

10.3969/j.issn.1674-0696.2015.06.04

2014-08-12；

2015-04-14

國家自然科學(xué)基金項(xiàng)目(51178396)

王 濤(1983—)， 男，四川南充人，博士研究生，主要從事橋梁結(jié)構(gòu)動(dòng)力學(xué)方面的研究。E-mail:7015294@qq.com。

U441;TB115

A

1674-0696(2015)06-019-08