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

王 濤，沈銳利

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

?

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

王 濤，沈銳利

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

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

橋梁工程；有限元程序；CR列式；幾何非線性；動力時程；斜拉索

0 引 言

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

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

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

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

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

1 CR列式方法基本原理

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

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

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

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

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

(1)

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

(2)

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

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

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

(3)

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

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

(4)

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

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

(5)

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

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

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

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

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

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

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

(6)

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

2 基于CR列式的非線性有限元動力時程積分

2.1 基本理論

筆者開發(fā)了基于CR列式的非線性Newmark-β有限元動力時程積分方法，其基本假定與普通的Newmark-β法相同，可以寫為：

(7)

(8)

其中各個參數(shù)具體表達式可以參見文獻[12]?？紤]幾何非線性時，結(jié)構(gòu)t+Δt時刻的振動方程為：

(9)

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

(10)

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

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

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

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

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

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

(11)

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

2.2 單元計算

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

(12)

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

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

(13)

式中：f為單元軸力。

(14)

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

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

(15)

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

(16)

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

(17)

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

(18)

2.3 算法流程

綜上所述，筆者基于CR列式的幾何非線性Newmark-β法有限元動力時程積分計算流程如圖4。

圖4 非線性有限元動力時程計算流程

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

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

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

3 算例驗證

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

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

3.1 算例1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 算例2

某實際斜拉橋拉索(圖10)的兩端點直線距離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個直桿單元。

圖10 斜拉索在端點位移激勵下的有限元模型

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

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

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

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

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

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

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

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

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

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

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

圖15 拉索各個時間點振動形狀(有放大)

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

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

4 結(jié) 論

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

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

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

[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é)報,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] 周凌遠,李喬.基于UL法的CR列式三維梁單元計算方法[J].西南交通大學(xué)學(xué)報,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è)計[M].北京：人民交通出版社,2006. Xu Rongqiao.Finite Element Method in Structural Analysis and MATLAB Programming [M].Beijing:China Communications Press,2006.

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

[14] 王濤,沈銳利.斜拉橋索-梁相關(guān)振動概念與研究方法初探[J].振動與沖擊,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é)基金項目(51178396)

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

U441;TB115

A

1674-0696(2015)06-019-08