王 濤，沈銳利

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

?

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

王 濤，沈銳利

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

基于CR列式非線性有限元計(jì)算理論，提出了針對(duì)桿、梁結(jié)構(gòu)的Newmark-β非線性有限元?jiǎng)恿r(shí)程積分算法,在動(dòng)力時(shí)程計(jì)算中，可以考慮結(jié)構(gòu)的幾何非線性效應(yīng)，通過(guò)計(jì)算單元伸長(zhǎng)并扣除結(jié)構(gòu)剛體位移得到單元的實(shí)際變形和內(nèi)力；闡述了非線性有限元?jiǎng)恿r(shí)程計(jì)算程序的計(jì)算原理,開(kāi)發(fā)了計(jì)算程序，建立有限元模型；通過(guò)算例與ANSYS的計(jì)算結(jié)果進(jìn)行了對(duì)比。結(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)。對(duì)幾何非線性，國(guó)內(nèi)外學(xué)者提出了多種方法。最常用的是以結(jié)構(gòu)變形前為參考，建立平衡方程的全拉格朗日法(TL列式法)和以結(jié)構(gòu)變形后為參考，建立平衡方程的更新拉格朗日法(UL列式法)。K.J.Bathe等[1-2]在大量研究基礎(chǔ)上，建立了三維梁?jiǎn)卧笪灰?、大轉(zhuǎn)動(dòng)、小應(yīng)變的UL列式和TL列式分析方法。陳政清等[3]改進(jìn)了K.J.Bathe建立的非線性梁?jiǎn)卧?，提高了UL列式計(jì)算效率。

對(duì)于大轉(zhuǎn)動(dòng)問(wèn)題，為了得到更為精確的結(jié)果，使用TL或UL列式往往需要增加荷載步，從而顯著增加了計(jì)算時(shí)間，且因涉及到大轉(zhuǎn)動(dòng)問(wèn)題，UL列式方法的推導(dǎo)過(guò)程往往都非常復(fù)雜，在實(shí)際編程的應(yīng)用中難度較大。國(guó)內(nèi)外廣大學(xué)者對(duì)該問(wèn)題其它解決途徑進(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]對(duì)不同單元幾何非線性CR列式提出了一致列式方法。K.Hsiao等[7]，C.A.Felippa等[8]對(duì)梁?jiǎn)卧腃R列式計(jì)算方法進(jìn)行了深入研究。周凌遠(yuǎn)等[9]將UL列式計(jì)算概念與CR列式方法結(jié)合起來(lái)，提出了基于UL法的CR列式三維梁?jiǎn)卧?jì)算方法。潘永仁[10]，唐茂林[11]基于CR列式幾何非線性計(jì)算方法開(kāi)發(fā)了有限元程序，運(yùn)用到大跨度懸索橋施工過(guò)程監(jiān)控，取得了良好的結(jié)果。

以上研究主要關(guān)注結(jié)構(gòu)靜力計(jì)算幾何非線性問(wèn)題。目前，在結(jié)構(gòu)動(dòng)力時(shí)程計(jì)算中，通常使用振型分解法或線性Newmark-β法來(lái)得到結(jié)構(gòu)的振動(dòng)響應(yīng)，對(duì)于結(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列式對(duì)單元內(nèi)力的計(jì)算方法，建立了考慮結(jié)構(gòu)幾何非線性的Newmark-β動(dòng)力時(shí)程計(jì)算方法，編制了計(jì)算程序，使用ANSYS驗(yàn)證了程序的正確性與可靠性。相對(duì)于一般基于UL列式幾何非線性計(jì)算方法的商業(yè)通用有限元軟件(如ANSYS等)，該方法計(jì)算速度快，力學(xué)概念簡(jiǎn)潔(通常不必計(jì)算、推導(dǎo)復(fù)雜的單元非線性大位移剛度矩陣)，適合于模擬柔性結(jié)構(gòu)(特別是索結(jié)構(gòu))的大幅度非線性振動(dòng)。

1 CR列式方法基本原理

以2維梁?jiǎn)卧獮槔?，在圖1(a)中給出了在結(jié)構(gòu)總體坐標(biāo)系XY中未變形的二維梁?jiǎn)卧?，單元受到總體坐標(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)的伸長(zhǎng)與轉(zhuǎn)角如圖1(b)。

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

(1)

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

(2)

對(duì)于梁?jiǎn)卧趩卧鲃?dòng)坐標(biāo)X1Y1中，單元節(jié)點(diǎn)的實(shí)際位移為單元伸長(zhǎng)后的長(zhǎng)度與初始長(zhǎng)度的差。單元節(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。

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

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

(3)

每一步迭代中梁?jiǎn)卧趩卧S動(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)可以看出，對(duì)于直桿單元，由于在流動(dòng)坐標(biāo)系中僅有X1軸方向的變形，使用CR列式方法計(jì)算得到的單元伸長(zhǎng)在幾何上是精確的。對(duì)于梁?jiǎn)卧?，得到的單元伸長(zhǎng)，通過(guò)扣除剛體轉(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可取為線彈性的桿、梁?jiǎn)卧仃嘯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)換矩陣的變化來(lái)反映的。由于CR列式為幾何精確方法，所以從理論上來(lái)講，荷載的分級(jí)加載的步數(shù)對(duì)最后結(jié)果無(wú)影響，差別來(lái)自于數(shù)值舍入的誤差。

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

(6)

對(duì)于式(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 基本理論

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

(7)

(8)

其中各個(gè)參數(shù)具體表達(dá)式可以參見(jiàn)文獻(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)的問(wèn)題。

式(10)左端中的總體節(jié)點(diǎn)內(nèi)力向量{R}是總體節(jié)點(diǎn)位移向量{u}的非線性函數(shù)。對(duì)于{R}可以按照?qǐng)D 2中的流程來(lái)計(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ù)，可以按照?qǐng)D3中的流程來(lái)計(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)換矩陣來(lái)更新[M]。對(duì)于阻尼，由于其復(fù)雜性，計(jì)算中不考慮[C]的變化。

可以使用牛頓迭代法來(lái)求解式(10)，當(dāng)?shù)词諗繒r(shí)，式(10)左右兩端是不相等的，計(jì)算兩端的差值就可以得到迭代的不平衡力。在迭代過(guò)程中使用的切線剛度矩陣為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)力剛度矩陣來(lái)組集。為了加快迭代計(jì)算收斂速度，式(11)也根據(jù)節(jié)點(diǎn)位置狀態(tài)來(lái)更新。

2.2 單元計(jì)算

筆者使用的梁?jiǎn)卧獜椥詣偠染仃嚍椋?/p>

(12)

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

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

(13)

式中：f為單元軸力。

(14)

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

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

(15)

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

(16)

桿、梁?jiǎn)卧淖鴺?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的位置與單元長(zhǎng)度l可計(jì)算得到：

(18)

2.3 算法流程

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

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

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

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

筆者只闡述了使用的關(guān)鍵技術(shù)方法，更多的關(guān)于有限元程序開(kāi)發(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)，開(kā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)畫(huà)圖形。

為了驗(yàn)證筆者程序的正確性與可靠性，將算例的結(jié)果與ANSYS計(jì)算結(jié)果進(jìn)行對(duì)比。程序計(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àn)的索-梁組合結(jié)構(gòu)有限元模型如圖5，主梁分為4個(gè)梁?jiǎn)卧鞴卜譃?0個(gè)直桿單元，彈性模量E、材料質(zhì)量密度ρ，單元截面積A，梁?jiǎn)卧箯潉偠菼z，桿單元初始軸力H，單元物理參數(shù)如表1。

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

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

考慮幾何非線性計(jì)算結(jié)構(gòu)在自重作用下的靜力構(gòu)型，得到ANSYS與筆者計(jì)算結(jié)果對(duì)比，如圖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í)間步長(zhǎng)為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)之間互相傳遞，呈“此消彼長(zhǎng)”的趨勢(shì)。

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

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

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

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é)果對(duì)比，如圖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í)間步長(zhǎng)取為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í)程，對(duì)比ANSYS與本文程序計(jì)算結(jié)果，如圖14。

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

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

筆者程序可以根據(jù)計(jì)算結(jié)果生成振動(dòng)時(shí)程的動(dòng)畫(huà)，可以更為直觀地觀察拉索非線性振動(dòng)的發(fā)展過(guò)程，但限于表達(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ì)地闡述了程序開(kāi)發(fā)過(guò)程。編制了有限元計(jì)算程序，通過(guò)算例的靜力與動(dòng)力計(jì)算與ANSYS結(jié)果的對(duì)比，驗(yàn)證了程序的正確、可靠性。

2)相對(duì)于ANSYS這樣使用UL列式計(jì)算結(jié)構(gòu)內(nèi)力的通用軟件，基于筆者算法開(kāi)發(fā)的有限元程序具有更加簡(jiǎn)潔、實(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)大撓度問(wèn)題內(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列式三維梁?jiǎn)卧?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] 唐茂林.大跨度懸索橋空間幾何非線性分析與軟件開(kāi)發(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

國(guó)家自然科學(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