王 濤， 劉德貴， 胡安杰

(西南科技大學(xué) 土木工程與建筑學(xué)院，四川 綿陽(yáng) 621000)

對(duì)于土木工程大跨度空間結(jié)構(gòu)，特別是大跨度橋梁，如：斜拉橋、懸索橋，考慮幾何非線性的有限元計(jì)算方法在結(jié)構(gòu)分析中的有著重要的應(yīng)用。

關(guān)于幾何非線性，最常用的是以結(jié)構(gòu)變形前為參考建立平衡方程的全量拉格朗日法(TL列式)和以結(jié)構(gòu)變形后為參考建立平衡方程的更新拉格朗日法(UL列式)。Bathe等[1]首先建立了三維梁?jiǎn)卧笪灰?、大轉(zhuǎn)動(dòng)、小應(yīng)變UL列式分析方法。陳政清等[2]詳細(xì)研究了基于UL列式的非線性有限元方法并將其運(yùn)用到了實(shí)際工程的計(jì)算中。

目前，大多數(shù)成熟的商業(yè)有限元軟件如ANSYS，通常都采用UL列式。但對(duì)于大轉(zhuǎn)動(dòng)問(wèn)題，為了在計(jì)算中得到更精確的結(jié)果，使用UL列式往往需要增加荷載分級(jí)加載步數(shù)，顯著增加了計(jì)算時(shí)間，且由于在非線性分析時(shí)材料本構(gòu)關(guān)系與單元的運(yùn)動(dòng)描述緊密耦合，UL列式方法的推導(dǎo)過(guò)程通常都非常繁瑣，在實(shí)際編程的應(yīng)用中難度相對(duì)較大，這也引起了國(guó)內(nèi)外廣大學(xué)者對(duì)該問(wèn)題其它解決途徑的研究興趣。算法簡(jiǎn)單可靠、編程方便、非線性計(jì)算收斂性與穩(wěn)定性好的流動(dòng)坐標(biāo)系方法(共旋坐標(biāo)系)，以下稱為CR列式法(Co-rotational Formulation)成為了結(jié)構(gòu)幾何非線性有限元計(jì)算方法研究的另一關(guān)注點(diǎn)。

Wempner[3]最早提出了基于CR列式幾何非線性的算法概念，即：認(rèn)為結(jié)構(gòu)是處于大變形、大位移、小應(yīng)變狀態(tài)，通過(guò)扣除剛體位移來(lái)得到有限元模型單元節(jié)點(diǎn)實(shí)際位移，該方法是一種“幾何精確”的算法，具有算法相對(duì)簡(jiǎn)單、力學(xué)概念明確等的優(yōu)點(diǎn)。Crisfield等[4]針對(duì)實(shí)體、殼、梁?jiǎn)卧膸缀畏蔷€性CR列式提出了一致列式方法。Felipppa等[5]對(duì)梁?jiǎn)卧腃R列式的計(jì)算方法進(jìn)行了深入研究。

Zhang等[6]使用CR列式非線性有限元方法計(jì)算了整體張拉結(jié)構(gòu)的大位移問(wèn)題，Liang等[7]基于CR列式有限元方法分析了結(jié)構(gòu)的非線性扭轉(zhuǎn)問(wèn)題，鄧?yán)^華等[8]開(kāi)發(fā)了基于CR列式并考慮徐變作用的非線性桿、梁?jiǎn)卧?，潘永仁等[9]基于CR列式幾何非線性計(jì)算方法開(kāi)發(fā)了有限元計(jì)算程序并運(yùn)用到了大跨度懸索橋施工過(guò)程監(jiān)控的計(jì)算中。

上述研究均主要關(guān)注的是CR列式在結(jié)構(gòu)靜力幾何非線性計(jì)算中的應(yīng)用。目前，在工程結(jié)構(gòu)的動(dòng)力時(shí)程計(jì)算中，很多時(shí)候使用振型分解法或線性的Newmark-β法來(lái)得到結(jié)構(gòu)的振動(dòng)響應(yīng)。對(duì)于結(jié)構(gòu)的小幅度振動(dòng)，線性計(jì)算方法尚能達(dá)到滿足工程要求的計(jì)算精度，但當(dāng)大跨度柔性結(jié)構(gòu)在3維空間發(fā)生較大幅度振動(dòng)時(shí)(特別是索結(jié)構(gòu))，振動(dòng)往往是非線性的，使用線性計(jì)算方法得到的計(jì)算結(jié)果很可能與實(shí)際情況偏差較大?？岛褴姷萚10]也在其研究中指出，對(duì)于像斜拉橋這樣的結(jié)構(gòu)，構(gòu)建高效、可靠的全橋非線性振動(dòng)計(jì)算大系統(tǒng)是橋梁工程學(xué)未來(lái)發(fā)展必須關(guān)注的研究方向。吳慶雄等[11]針對(duì)索梁結(jié)構(gòu)非線性振動(dòng)建立了有限元計(jì)算方法，開(kāi)發(fā)了計(jì)算程序，但文中算法是建立在結(jié)構(gòu)靜力平衡狀態(tài)上的，計(jì)算中不能直接考慮結(jié)構(gòu)非線性振動(dòng)時(shí)的重力變化，也沒(méi)有考慮端點(diǎn)強(qiáng)制位移激勵(lì)作用下的非線性振動(dòng)計(jì)算方法。曹九發(fā)等[12]建立了使用完全拉格朗日方法(TL列式)的非線性動(dòng)力有限元程序，分析了大型風(fēng)力機(jī)的幾何非線性動(dòng)態(tài)響應(yīng)，但TL列式始終以結(jié)構(gòu)的初始構(gòu)型為參考，很多情況下，為了簡(jiǎn)化推到過(guò)程，計(jì)算使用的非線性切線剛度矩陣本身也是近似的，當(dāng)柔性結(jié)構(gòu)具有較大的振動(dòng)位移時(shí)，計(jì)算結(jié)果可能與實(shí)際情況呈現(xiàn)較大偏差。

為了在動(dòng)力時(shí)程計(jì)算中模擬工程結(jié)構(gòu)發(fā)生大幅振動(dòng)時(shí)的幾何非線性效應(yīng)，本文基于CR列式首先建立了考慮初始應(yīng)力的3 維空間中非線性桿、梁?jiǎn)卧缓箝_(kāi)發(fā)了考慮幾何非線性的CR列式Newmark-β動(dòng)力時(shí)程計(jì)算方法，編制了通用的非線性動(dòng)力計(jì)算程序(使用桿、梁?jiǎn)卧?，通過(guò)與ANSYS計(jì)算對(duì)比，驗(yàn)證了本文算法具有較快的計(jì)算速度且編制的程序具有良好的可靠性、穩(wěn)定性，適用于作為大跨度橋梁結(jié)構(gòu)(斜拉橋、懸索橋)幾何非線性振動(dòng)研究的基礎(chǔ)計(jì)算方法。

2 基于CR列式桿梁結(jié)構(gòu)幾何非線性計(jì)算原理

2.1 三種坐標(biāo)系統(tǒng)

根據(jù)潘永仁等人的論述，本文認(rèn)為，基于CR列式非線性有限元算法的核心問(wèn)題在于建立與有限元模型單元相關(guān)的3種坐標(biāo)系統(tǒng)(即：結(jié)構(gòu)建模總體坐標(biāo)系、單元隨動(dòng)坐標(biāo)系、單元端截面隨動(dòng)坐標(biāo)系)以及處理它們之間的關(guān)系。由于直桿單元屬于梁?jiǎn)卧耐嘶闆r，所以這里以3維梁?jiǎn)卧獮槔齺?lái)闡述桿、梁?jiǎn)卧?種坐標(biāo)系之間的關(guān)系。

(1)結(jié)構(gòu)總體坐標(biāo)系

如圖1所示，結(jié)構(gòu)總體坐標(biāo)系用xi(i=1, 2, 3)表示(對(duì)應(yīng)X,Y,Z軸)，其基矢量為ei(i=1, 2, 3), 總體坐標(biāo)系始終不變，結(jié)構(gòu)的建模與最終計(jì)算得到的模型節(jié)點(diǎn)位移，均依據(jù)總體坐標(biāo)系來(lái)描述。

(2)單元隨動(dòng)坐標(biāo)系

如圖1所示，為了得到空間梁?jiǎn)卧膯卧植孔鴺?biāo)與整體坐標(biāo)系的關(guān)系。這里定義梁?jiǎn)卧膯卧S動(dòng)坐標(biāo)系。單元隨動(dòng)坐標(biāo)系描述了整個(gè)梁?jiǎn)卧?維空間中的位置。

圖1 總體坐標(biāo)系與梁?jiǎn)卧S動(dòng)坐標(biāo)系Fig.1 The global coordinate system and co-rotationalcoordinate system of beam element

(3)單元端截面隨動(dòng)坐標(biāo)系

圖2 總體坐標(biāo)系與梁?jiǎn)卧私孛骐S動(dòng)坐標(biāo)系Fig.2 The global coordinate system and co-rotationalcoordinate systems of the cross section of beam element

1.2 單元變形后隨動(dòng)坐標(biāo)系的計(jì)算

對(duì)于大變形梁?jiǎn)卧?，要得到正確的單元節(jié)點(diǎn)力，關(guān)鍵在于計(jì)算單元在3維空間中的伸長(zhǎng)，以及梁端截面的轉(zhuǎn)角和單元變形后的單元隨動(dòng)坐標(biāo)系。

(1)歐拉有限轉(zhuǎn)動(dòng)公式

對(duì)于3維空間矢量旋轉(zhuǎn)，本文使用歐拉有限轉(zhuǎn)動(dòng)公式[9]。如圖3所示，矢量R的初始位置用OP來(lái)表示，按照右手定則，沿著轉(zhuǎn)動(dòng)軸ON逆時(shí)針旋轉(zhuǎn)了θ角成為R′，ON方向用單位矢量n來(lái)表示，垂直于轉(zhuǎn)動(dòng)軸過(guò)N點(diǎn)的平面如圖3(b)所示。

圖3 有限轉(zhuǎn)動(dòng)矢量示意圖Fig.3 The schematic figure of finite rotation vector

根據(jù)圖3中各個(gè)矢量的關(guān)系可以求出R′：

ON=n(n·R)
NP=R-n(n·R)
R′=ON+NV+VQ

即：

R′=Rcosθ+n(n·R)(1-cosθ)+(n×R)sinθ

(1)

(2)單元端截面隨動(dòng)坐標(biāo)系的計(jì)算

已知在t至t+Δt時(shí)刻，在結(jié)構(gòu)坐標(biāo)系中單元節(jié)點(diǎn)發(fā)生了轉(zhuǎn)動(dòng)位移矢量Δθji，平動(dòng)位移矢量ΔUji，j=1,2;i=1, 2, 3。

根據(jù)Δθji由數(shù)學(xué)定義可知轉(zhuǎn)動(dòng)角度值為：

(2)

則各節(jié)點(diǎn)的轉(zhuǎn)動(dòng)軸矢量方向?yàn)椋?/p>

(3)

(3)單元隨動(dòng)坐標(biāo)系的計(jì)算

設(shè)t時(shí)刻在總體坐標(biāo)系下，單元節(jié)點(diǎn)坐標(biāo)為tXji，則可以得到t+Δt時(shí)刻節(jié)點(diǎn)坐標(biāo)為：

t+ΔtXji=tXji+ΔUji,j=1, 2;i=1, 2, 3

(4)

(5)

(6)

1.3 單元隨動(dòng)坐標(biāo)系中單元變形的計(jì)算

在t+Δt時(shí)刻單元變形可以分解為：彎曲變形，扭轉(zhuǎn)變形，伸長(zhǎng)變形。

(1)彎曲變形

(7)

(2)扭轉(zhuǎn)變形

(3)伸長(zhǎng)變形

在潘永仁的研究中，將單元的伸長(zhǎng)變形定義為薄膜變形，通過(guò)定義單元變形曲線并計(jì)算曲線弧長(zhǎng)變化來(lái)得到單元的薄膜變形。

本文認(rèn)為，對(duì)于一般的幾何非線性計(jì)算，根據(jù)小應(yīng)變假設(shè)，從直觀上來(lái)看可以由節(jié)點(diǎn)的位置直接計(jì)算得到單元的伸長(zhǎng)變形值。

已知t時(shí)刻單元節(jié)點(diǎn)位置為tXji，則可以計(jì)算得到t時(shí)刻單元長(zhǎng)度為

(8)

根據(jù)式(4)得到總體坐標(biāo)系下t+Δt時(shí)刻節(jié)點(diǎn)坐標(biāo)t+ΔtXji，與式(8)方法相同計(jì)算t+Δt時(shí)刻單元長(zhǎng)度t+Δtl可以得到

Δl=t+Δtl-tl

(9)

Δl即為t+Δt時(shí)刻單元的伸長(zhǎng)值。

1.4 節(jié)點(diǎn)位移計(jì)算

(1)單元節(jié)點(diǎn)力的計(jì)算

根據(jù)1.2節(jié)與1.3節(jié)中的計(jì)算原理，可以得到在單元隨動(dòng)坐標(biāo)系中單元變形后t+Δt時(shí)刻的單元節(jié)點(diǎn)位移向量ue，由于單元隨動(dòng)坐標(biāo)系原點(diǎn)在單元節(jié)點(diǎn)1上，所以u(píng)e可寫為

ue=[0 0 0θ11θ12θ13Δl0 0θ21θ22θ23]

(10)

單元隨動(dòng)坐標(biāo)系中單元實(shí)際內(nèi)力為：

fe=keue

(11)

式中：fe為單元內(nèi)力的節(jié)點(diǎn)力向量；ke為單元彈性剛度矩陣。對(duì)于直桿單元，流動(dòng)坐標(biāo)系內(nèi)單元的位移只有沿軸向的伸長(zhǎng)而無(wú)轉(zhuǎn)動(dòng)，在流動(dòng)坐標(biāo)系ke內(nèi)取為桿單元的線彈性剛度矩陣即為“幾何精確的”； 對(duì)于梁?jiǎn)卧?，一般情況下考慮結(jié)構(gòu)幾何非線性計(jì)算時(shí)，通過(guò)CR列式方法，扣除單元?jiǎng)傮w轉(zhuǎn)動(dòng)后，可以認(rèn)為梁?jiǎn)卧翘幱诖笞冃?、大轉(zhuǎn)動(dòng)、小應(yīng)變狀態(tài)的，那么在單元的局部流動(dòng)坐標(biāo)系內(nèi)，則可以認(rèn)為力與位移的關(guān)系是近似線性的，ke可以近似取為梁?jiǎn)卧木€彈性剛度矩陣。

(2)單元坐標(biāo)轉(zhuǎn)換矩陣與應(yīng)力剛度矩陣

當(dāng)單元中存在初始軸力時(shí)，單元的轉(zhuǎn)向與側(cè)向變形會(huì)導(dǎo)致單元產(chǎn)生應(yīng)力剛度矩陣，在桿梁?jiǎn)卧木€性計(jì)算中，單元應(yīng)力剛度矩陣通常是必須考慮的，而在非線性計(jì)算中，單元應(yīng)力剛度矩陣僅用來(lái)組集總體切線剛度矩陣，加快迭代收斂。3維空間桿、梁?jiǎn)卧獞?yīng)力剛度矩陣的表達(dá)式可參考文獻(xiàn)[13-14]。

(3)節(jié)點(diǎn)位移計(jì)算

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

圖4 有限元模型總體節(jié)點(diǎn)內(nèi)力向量計(jì)算流程Fig.4 The calculation process of the global nodeinternal force vector of FEM modal

結(jié)構(gòu)產(chǎn)生大變形時(shí)，有限元模型的總體節(jié)點(diǎn)外力向量可表示為F(u)，也為位移的非線性函數(shù)。直接作用在節(jié)點(diǎn)上的外力在通常情況下可以認(rèn)為其不變，如果發(fā)生改變也可以使用分級(jí)加載實(shí)現(xiàn)外力的變化。而重力加速度、單元初始應(yīng)力的造成的等效外力，必然會(huì)隨著結(jié)構(gòu)的大變形而改變，因此本文計(jì)算中考慮了重力與單元初始應(yīng)力的非線性時(shí)變效應(yīng)?？傮w節(jié)點(diǎn)外力向量計(jì)算流程，如圖5所示。

圖5 有限元模型總體節(jié)點(diǎn)外力向量計(jì)算流程Fig.5 The calculation process of the global nodeexternal force vector of FEM modal

根據(jù)圖4與圖5的計(jì)算流程，可以看出，在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)，有如下關(guān)系：

R(u)=F(u)

(12)

靜力計(jì)算的最終目的就是得到結(jié)構(gòu)在外力作用下達(dá)到靜力平衡狀態(tài)時(shí)總體坐標(biāo)系下有限元模型各個(gè)節(jié)點(diǎn)的位移向量u。對(duì)于式(12)依據(jù)CR列式計(jì)算原理，一般使用Newton-Raphson迭代法計(jì)算，迭代中所需要的總體切線剛度矩陣可以通過(guò)疊加當(dāng)前時(shí)間步中結(jié)構(gòu)總體彈性剛度矩陣與應(yīng)力剛度矩陣來(lái)得到。

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

2.1 基本理論

在3維桿、梁?jiǎn)卧蔷€性靜力計(jì)算的基礎(chǔ)上，本文開(kāi)發(fā)了基于CR列式的非線性Newmark-β有限元?jiǎng)恿r(shí)程積分方法，在算法中也考慮了結(jié)構(gòu)初始應(yīng)力的影響。其基本假定與普通的Newmark-β法相同：

(13a)

(13b)

(14)

將式(13)代入式(14)可以得到：

(15)

這樣，根據(jù)Newmark-β的基本原理，非線性有限元?jiǎng)恿r(shí)程積分就被轉(zhuǎn)化為在每一個(gè)時(shí)間步內(nèi)求解非線性方程組(15)的問(wèn)題，可以采用CR列式方法原理，使用Newton-Raphson法來(lái)進(jìn)行平衡迭代計(jì)算。

式(15)左端中的總體節(jié)點(diǎn)內(nèi)力向量R是總體節(jié)點(diǎn)位移向量u的非線性函數(shù)。對(duì)于R可以按照?qǐng)D 4中的流程來(lái)計(jì)算 。

在考慮結(jié)構(gòu)幾何非線性時(shí)，式(15)右端中，F(xiàn)表示結(jié)構(gòu)承受的總體節(jié)點(diǎn)外力向量(如：重力、結(jié)構(gòu)上施加的外力)，F(xiàn)e表示結(jié)構(gòu)初始單元應(yīng)力造成的等效總體節(jié)點(diǎn)外力向量(與ANSYS中初應(yīng)變導(dǎo)致單元等效外力的概念相同，在本文有限元程序建模中可以直接定義單元初始軸力，拉力為正，壓力為負(fù))。F，F(xiàn)e都是位移u的非線性函數(shù)，可以按照?qǐng)D5中的流程來(lái)計(jì)算。

式(15)左右兩端還存在與總體質(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。由于阻尼的復(fù)雜性且在一般情況下結(jié)構(gòu)阻尼的變化對(duì)振動(dòng)的影響不大，在本文計(jì)算中不考慮C的改變。

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

K=a0M+a1C+KK

(16)

式中：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ì)算收斂速度，式(16)也根據(jù)節(jié)點(diǎn)位置狀態(tài)來(lái)更新。

2.2 算法流程

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

圖6 非線性有限元?jiǎng)恿r(shí)程計(jì)算流程圖Fig.6 Calculation flow graph of nonlinear time-history FEM

3 算例驗(yàn)證

依據(jù)前文所述原理，使用MATLAB數(shù)值計(jì)算平臺(tái)，本文開(kāi)發(fā)了基于CR列式的3維空間結(jié)構(gòu)非線性動(dòng)力有限元程序，程序主要實(shí)現(xiàn)的功能如下：①可以考慮幾何非線性計(jì)算結(jié)構(gòu)靜力狀態(tài)；②可以在得到結(jié)構(gòu)發(fā)生非線性大變形的靜力狀態(tài)后，計(jì)算結(jié)構(gòu)的動(dòng)力特性；③可以在計(jì)算中考慮初始應(yīng)力與自重的時(shí)變效應(yīng)，計(jì)算結(jié)構(gòu)在外激勵(lì)或強(qiáng)制位移激勵(lì)作用下的非線性動(dòng)力時(shí)程響應(yīng)；④可以根據(jù)非線性動(dòng)力時(shí)程計(jì)算結(jié)果繪制輸出結(jié)構(gòu)非線性振動(dòng)的動(dòng)畫。

3.1 算例1 彎曲梁3維空間靜力計(jì)算

本文中CR列式3維非線性動(dòng)力時(shí)程計(jì)算是以非線性靜力計(jì)算為基礎(chǔ)的，所以，這里首先使用一個(gè)經(jīng)典的靜力模型來(lái)驗(yàn)證算法的正確性。

如圖7所示，45°彎曲梁空間彎扭幾何非線性大變形計(jì)算。與參考文獻(xiàn)[9]中的結(jié)構(gòu)相同，在XYZ坐標(biāo)系中，梁的B端固結(jié)位于坐標(biāo)系原點(diǎn)，A端不受約束且受到沿Y方向的力，梁共分為8個(gè)等長(zhǎng)度的空間梁?jiǎn)卧?，各個(gè)節(jié)點(diǎn)坐標(biāo)按照?qǐng)A弧曲線計(jì)算。材料彈性模型為E=107，泊松比u=0，剪切模量G=E/2(1+u) ,圓弧半徑R=100，截面尺寸為的矩形1.0×1.0。端點(diǎn)力P在計(jì)算迭代過(guò)程中始終保持沿Y方向。非線性靜力計(jì)算中本文程序與ANSYS均設(shè)置12級(jí)分級(jí)加載。

圖7 彎曲梁3維空間模型Fig.7 The mode of a curving beam in 3D space

圖8中給出了P=300,P=600時(shí)，結(jié)構(gòu)在靜力作用下發(fā)生大幅度變形后的構(gòu)型。

圖8 彎曲梁3維空間變形Fig.8 The displacement of the curving beam in 3D space

已知彎曲梁在變形前，梁端節(jié)點(diǎn)A的位置坐標(biāo)為x=70.7,y=0,z=29.3。使用本文程序計(jì)算得到結(jié)構(gòu)變形后的坐標(biāo)位置與參考文獻(xiàn)[9]和ANSYS計(jì)算得到節(jié)點(diǎn)A的位置結(jié)果對(duì)比，如表1所示。

表1 計(jì)算結(jié)果對(duì)比

由表1可以看出，由于都采用基于CR列式的幾何非線性靜力計(jì)算方法，本文計(jì)算結(jié)果與參考文獻(xiàn)差別很小，本文認(rèn)為這是由于程序設(shè)計(jì)的細(xì)節(jié)差別造成的。由于結(jié)構(gòu)的變形很大，且ANSYS中幾何非線性計(jì)算使用的是UL列式，所以ANSYS結(jié)果與本文以及文獻(xiàn)[9]的差別相對(duì)較明顯一些。

3.2 算例2：索-梁組合結(jié)構(gòu)3維空間非線性振動(dòng)計(jì)算

橋梁工程中常見(jiàn)的索-梁組合結(jié)構(gòu)簡(jiǎn)化模型[15]，有限元模型的幾何尺寸，如圖9所示。

圖9 索-梁組合結(jié)構(gòu)有限元模型Fig.9 Cable-Beam structure FEM modal

總體坐標(biāo)系為XYZ，建立拉索模型局部坐標(biāo)系為X1Y1Z1，Z1與Z軸指向相同。梁分為4個(gè)3維梁?jiǎn)卧?，拉索共分?0個(gè)3維直桿單元。不計(jì)泊松比，結(jié)構(gòu)單元的各個(gè)物理參數(shù)為：彈性模量E, Pa、剪切模量G, Pa、材料質(zhì)量密度ρ, kg/m3、單元截面積A, m2、梁?jiǎn)卧箯潙T性矩Iz, m4,Iy, m4、抗扭慣性矩Ix, m4、桿單元初始軸力H, N，如表2所示。

表2 結(jié)構(gòu)單元物理參數(shù)

為了驗(yàn)證本文程序的正確性與可靠性，通過(guò)算例的結(jié)果與ANSYS對(duì)比。ANSYS中使用Beam4梁?jiǎn)卧M主梁，使用Link8桿單元模擬拉索，這兩種單元的都為3維單元[16]。本文程序與ANSYS計(jì)算均使用一致質(zhì)量矩陣。

程序計(jì)算中的重力加速度均取為G=9.8 m/s2。計(jì)算中不設(shè)置結(jié)構(gòu)阻尼。由文獻(xiàn)[16]可知ANSYS計(jì)算中默認(rèn)的算法阻尼不為零，即：Newmark-β法中參數(shù)γ=0.505，在本文程序計(jì)算中使用相同的設(shè)置。

首先，使用本文有限元程序考慮幾何非線性計(jì)算結(jié)構(gòu)在自重作用下的靜力構(gòu)型，然后計(jì)算動(dòng)力特性，在總體坐標(biāo)系XYZ下得到前6階模態(tài)計(jì)算結(jié)果，如圖10所示。

圖10 結(jié)構(gòu)前6階模態(tài)Fig.10 The first six modes of the structure

非線性靜力計(jì)算結(jié)果得到結(jié)構(gòu)在自重下拉索的平均軸力(各單元軸力之和除以單元數(shù)量)為50.65 kN，根據(jù)文獻(xiàn)[15]中的計(jì)算公式，單獨(dú)計(jì)算拉索的自振頻率，可以得到拉索在X1Y1平面內(nèi)1階自振頻率為3.166 5 Hz，拉索在X1Z1平面內(nèi)的1階自振頻率為3.153 9 Hz。

觀察圖10 (第2階振型)可以看出，由于拉索局部的1階自振頻率接近整體結(jié)構(gòu)的1階自振頻率，根據(jù)非線性振動(dòng)理論， 整體結(jié)構(gòu)發(fā)生1階振動(dòng)時(shí)，在端點(diǎn)位移激勵(lì)(節(jié)點(diǎn)5)沿拉索局部坐標(biāo)系Y1方向的分量作用下，拉索會(huì)發(fā)生1階非線性1∶1主共振。

然后，在節(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)在平面內(nèi)有自重狀態(tài)下的振動(dòng)時(shí)程，如圖11所示。

圖11 索-梁組合結(jié)構(gòu)面內(nèi)非線性振動(dòng)ANSYS與本文程序計(jì)算結(jié)果Fig.11 The results of the Cable-Beam structure in-plane nonlinear vibration which calculated by ANSYS and the program of this thesis

從圖11可看出，在釋放節(jié)點(diǎn)力后結(jié)構(gòu)發(fā)生了振動(dòng)，拉索在梁的帶動(dòng)下發(fā)生了1∶1主共振，拉索端部(節(jié)點(diǎn)5)沿整體坐標(biāo)系Y方向0.01 m幅值的位移激勵(lì)，導(dǎo)致拉索中點(diǎn)(節(jié)點(diǎn)11)沿局部坐標(biāo)系Y1方向發(fā)生了0.06 m的振動(dòng)。結(jié)構(gòu)振動(dòng)體現(xiàn)了“拍振[17]”的非線性振動(dòng)性質(zhì)，振動(dòng)能量在拉索與整體結(jié)構(gòu)之間互相傳遞，呈現(xiàn)“此消彼長(zhǎng)”的趨勢(shì)。符合非線性振動(dòng)的理論預(yù)期。

如圖11所示，由于非線性動(dòng)力時(shí)程計(jì)算中考慮了重力的直接作用，本文與文獻(xiàn)[11]不同，本文計(jì)算結(jié)果中振動(dòng)位移初始數(shù)值不為零，結(jié)構(gòu)模型在自重作用下非線性靜力計(jì)算后的初始位置(偏離y=0位置)為振動(dòng)平衡位置。有限元模型主梁端部節(jié)點(diǎn)5在自重作用下的靜力位移較小(約0.001 m)，偏離較不明顯，而拉索中點(diǎn)節(jié)點(diǎn)11(模型自重初始靜力位移約為0.02 m)的振動(dòng)時(shí)程圖可清晰地看出這個(gè)計(jì)算結(jié)果。

在單元節(jié)點(diǎn)5上施加簡(jiǎn)諧力外激勵(lì)P=P0sin(ωt)，其中P0取為1.0 kN，ω對(duì)應(yīng)的頻率取為3.16 Hz。由于簡(jiǎn)諧力與結(jié)構(gòu)面內(nèi)的1階自振頻率接近，所以，在外激勵(lì)作用下整體結(jié)構(gòu)會(huì)發(fā)生共振，從而導(dǎo)致拉索發(fā)生1∶1主共振。得到拉索1/2點(diǎn)的沿局部坐標(biāo)系Y1方向振動(dòng)時(shí)程計(jì)算結(jié)果，如圖12所示。

圖12 索-梁組合結(jié)構(gòu)在外激勵(lì)作用下拉索1/2點(diǎn)的沿Y1方向振動(dòng)時(shí)程圖、頻譜圖Fig.12 The time-history curves and spectrogram of the 1/2 point of the cable nonlinear vibrationin Cable-Beam structure in Y1direction under external excitation

從圖12中可以看出，結(jié)構(gòu)在簡(jiǎn)諧激勵(lì)下發(fā)生了共振，拉索振幅迅速增加(最大約0.3 m)但由于幾何非線性的作用，拉索振幅不會(huì)持續(xù)增加，根據(jù)文獻(xiàn)[17]的理論描述，這是由“振動(dòng)硬化”現(xiàn)象導(dǎo)致的，即：幾何非線性造成結(jié)構(gòu)的振動(dòng)頻率隨著振幅增大而增加。從計(jì)算結(jié)果頻譜圖中可以看出拉索的振動(dòng)包含了比自振頻率更高頻的成分。同事，計(jì)算結(jié)果時(shí)程圖中觀察到了非線性振動(dòng)造成拉索振幅隨時(shí)間“漲落”的更為明顯的“拍振”現(xiàn)象。

這里如果使用線性動(dòng)力時(shí)程計(jì)算，由于結(jié)構(gòu)剛度不變，在結(jié)構(gòu)發(fā)生共振時(shí)，理論上振幅會(huì)持續(xù)增加，直致小阻尼限制的上限，會(huì)遠(yuǎn)大于非線性振動(dòng)幅值上限，也不會(huì)發(fā)生非線性振動(dòng)特有的振動(dòng)硬化現(xiàn)象。所以，柔性結(jié)構(gòu)的大幅振動(dòng)必須考慮幾何非線性效應(yīng)。

對(duì)比圖11與12中使用ANSYS與本文程序計(jì)算得到的振動(dòng)時(shí)程圖，二者在振動(dòng)發(fā)展初期幾乎是一致的，由于ANSYS非線性有限元計(jì)算使用的是更新拉格朗日格式(UL列式)與本文CR列式有一定的區(qū)別，所以，兩者的計(jì)算結(jié)果存在細(xì)微差別，隨著時(shí)間增加，動(dòng)力時(shí)程積分的累加作用會(huì)導(dǎo)致差別相對(duì)小幅度增加。二者的時(shí)程計(jì)算結(jié)果都符合非線性振動(dòng)的理論預(yù)期。因此，通過(guò)與ANSYS的對(duì)比驗(yàn)證，筆者認(rèn)為本文程序的計(jì)算結(jié)果是正確、可靠的。

為了討論索-梁組合結(jié)構(gòu)的面外非線性振動(dòng)，在圖9模型節(jié)點(diǎn)5上，施加沿Z方向的橫向力P=25 kN，靜力計(jì)算后釋放節(jié)點(diǎn)力，動(dòng)力時(shí)程積分取時(shí)間步長(zhǎng)為0.01 s，使用本文程序計(jì)算結(jié)構(gòu)在有自重狀態(tài)下的振動(dòng)時(shí)程，得到梁上節(jié)點(diǎn)5與拉索1/2點(diǎn)(節(jié)點(diǎn)11)沿面外Z方向的振動(dòng)時(shí)程，如圖13所示。

圖13 索-梁組合結(jié)構(gòu)沿面外Z方向非線性振動(dòng)時(shí)程圖、頻譜圖Fig.13 The time-history curves and spectrograms of the cable-beam structure nonlinear vibration in out-plane Z direction

圖14 索-梁組合結(jié)構(gòu)非線性振動(dòng)時(shí)程狀態(tài)Fig.14 The time history states of nonlinearvibration of the cable-beam structure

同樣觀察圖10(第1階振型)可以知道，由于拉索的面外1階自振頻率與整體結(jié)構(gòu)接近，拉索在梁的帶動(dòng)下發(fā)生了面外主共振(沿Z方向)拉索發(fā)生了相對(duì)較大幅度的振動(dòng)。梁的面外抗彎剛度較大，所以“拍振”現(xiàn)象對(duì)梁的振幅漲落相對(duì)較不明顯。由于振動(dòng)硬化現(xiàn)象的作用，結(jié)構(gòu)響應(yīng)頻率中(4.0 Hz)存在大于面外自振頻率(3.148 Hz)的成分。

使用MATLAB可視化動(dòng)畫技術(shù)，我們可以更加直觀地觀察到結(jié)構(gòu)在3維空間的非線性振動(dòng)狀態(tài)，但限于表達(dá)方式，本文這里僅根據(jù)計(jì)算結(jié)果列出結(jié)構(gòu)發(fā)生面外振動(dòng)(橫向)在0～1.2 s內(nèi)的非線性振動(dòng)時(shí)程狀態(tài)(振動(dòng)位移放大5倍)，如圖14所示。

從圖14中可以看出，由于非線性振動(dòng)效應(yīng)，計(jì)算結(jié)果中觀察到了拉索的振動(dòng)滯后現(xiàn)象，即：梁向左端振動(dòng)、拉索向右端振動(dòng)的狀態(tài)(如第0.65 s)。這是更加符合實(shí)際情況的計(jì)算結(jié)果。

3.3 算例3:拉索3維空間回旋運(yùn)動(dòng)計(jì)算

如圖 15所示，一根緊繃的水平布置拉索，AB兩端固定，受到重力作用。在拉索端點(diǎn)A上作用圓弧型的位移激勵(lì)，拉索發(fā)生空間運(yùn)動(dòng)。

圖15 拉索受端點(diǎn)位移激勵(lì)示意圖Fig.15 The schematic figure of the cablewith displacement stimulation at the end point

已知拉索兩端的長(zhǎng)度為l=100 m(建模長(zhǎng)度)，彈性模量E=2.01×1011Pa, 截面積A=0.08 m2，拉索的初始軸力為H=4 000 kN，重力加速度G=9.8 m/s2，拉索分為20個(gè)3維直桿單元。

設(shè)拉索AB端水平固定，考慮幾何非線性，靜力計(jì)算得到拉索的自重狀態(tài)，垂度為-0.194 7 m，然后，計(jì)算自振頻率得到拉索面內(nèi)(XY平面)為1.257 7 Hz 面外自振頻率(XZ平面)為1.252 8 Hz。

在得到拉索的自重狀態(tài)后，在拉索端點(diǎn)施加圓弧軌跡位移激勵(lì)，如圖15所示，設(shè)端點(diǎn)A圓弧運(yùn)動(dòng)半徑R=0.005 m，在ZY坐標(biāo)系中Z方向分量為z=Rsin(ωt)，Y方向運(yùn)動(dòng)分量為Rcos(ωt)，設(shè)置拉索端點(diǎn)A圓弧運(yùn)動(dòng)頻率為1.253 Hz，則ω=1.253×2×π=7.872 8 rad/s。

使用強(qiáng)制位移施加端點(diǎn)位移激勵(lì)，取動(dòng)力時(shí)程積分步長(zhǎng)為0.02 s，計(jì)算2 000步。為了使“拍振”減小，振動(dòng)盡快進(jìn)入穩(wěn)定狀態(tài)，須設(shè)置較為明顯的阻尼作用，不設(shè)置結(jié)構(gòu)阻尼，設(shè)置算法阻尼γ=0.55。得到拉索中點(diǎn)在3維空間振動(dòng)的運(yùn)動(dòng)軌跡曲線，如圖16所示。

圖16 拉索中點(diǎn)3維空間非線性振動(dòng)曲線Fig.16 The nonlinear vibration curveof the 1/2 point of the cable in 3D space

從圖16中可以看出，微小的端點(diǎn)位移激勵(lì)，導(dǎo)致拉索發(fā)生較大幅度的空間振動(dòng)，拉索中點(diǎn)的運(yùn)動(dòng)軌跡呈明顯的圓曲線狀態(tài)。

造成上述現(xiàn)象的原因?yàn)椋河捎诶鞯乃髁^大，在自重作用下，拉索面外1階自振頻率仍然接近面內(nèi)1階自振頻率，而拉索端點(diǎn)位移激勵(lì)沿圓周運(yùn)動(dòng)且頻率接近拉索的面內(nèi)與面外1階自振頻率，因此，在端點(diǎn)位移激勵(lì)下拉索發(fā)生了1階面內(nèi)-面外耦合的1∶1主共振，導(dǎo)致了較大幅度的3空間回旋運(yùn)動(dòng)，即：柔性索的“跳繩”運(yùn)動(dòng)現(xiàn)象。

4 結(jié) 語(yǔ)

(1)開(kāi)發(fā)了3維空間桿梁?jiǎn)卧贑R列式的非線性Newmark-β動(dòng)力時(shí)程有限元算法，闡述了程序計(jì)算原理，編制有限元計(jì)算程序，通過(guò)算例與ANSYS計(jì)算結(jié)果對(duì)比，驗(yàn)證了程序的正確性與可靠性。

(2)基于本文算法開(kāi)發(fā)的有限元程序具有通用性，可以計(jì)算桿梁結(jié)構(gòu)在3維空間中外力作用下的非線性振動(dòng)，也可以計(jì)算結(jié)構(gòu)在強(qiáng)制位移激勵(lì)下的非線性振動(dòng)，可以在動(dòng)力時(shí)程計(jì)算中考慮初始應(yīng)力與重力隨結(jié)構(gòu)大幅振動(dòng)的時(shí)變效應(yīng)。

(3)經(jīng)過(guò)算例驗(yàn)證表明，本文的CR列式非線性動(dòng)力時(shí)程算法具有簡(jiǎn)潔、可靠、高效的特點(diǎn)。在筆者計(jì)算機(jī)上，采用同樣的計(jì)算設(shè)置，使用ANSYS進(jìn)行非線性動(dòng)力時(shí)程計(jì)算，算例1耗時(shí)約為3s，算例2約為150 s，算例3約為120 s，而使用本文程序耗時(shí)分別約為1 s、30 s、23 s。

(4)非線性有限元?jiǎng)恿r(shí)程計(jì)算能更好地反映結(jié)構(gòu)的受力細(xì)節(jié)，為柔性結(jié)構(gòu)的非線性振動(dòng)研究提供了較為完善的解決方案。對(duì)于橋梁工程專業(yè)，自主開(kāi)發(fā)計(jì)算程序更有利于整合本專業(yè)的各個(gè)算法，如：風(fēng)-車-橋耦合振動(dòng)算法。計(jì)劃在后續(xù)的研究中完善本程序，優(yōu)化計(jì)算效率，進(jìn)一步提高非線性動(dòng)力時(shí)程積分的計(jì)算速度并將其運(yùn)用到實(shí)際大跨度斜拉橋、懸索橋3維空間全橋模型非線性振動(dòng)的計(jì)算分析中。