王 濤， 劉德貴， 胡安杰

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

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

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

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

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

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

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

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

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

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

根據(jù)潘永仁等人的論述，本文認(rèn)為，基于CR列式非線性有限元算法的核心問題在于建立與有限元模型單元相關(guān)的3種坐標(biāo)系統(tǒng)(即：結(jié)構(gòu)建?？傮w坐標(biāo)系、單元隨動坐標(biāo)系、單元端截面隨動坐標(biāo)系)以及處理它們之間的關(guān)系。由于直桿單元屬于梁單元的退化情況，所以這里以3維梁單元為例來闡述桿、梁單元中3種坐標(biāo)系之間的關(guān)系。

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

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

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

如圖1所示，為了得到空間梁單元的單元局部坐標(biāo)與整體坐標(biāo)系的關(guān)系。這里定義梁單元的單元隨動坐標(biāo)系。單元隨動坐標(biāo)系描述了整個梁單元在3維空間中的位置。

圖1 總體坐標(biāo)系與梁單元隨動坐標(biāo)系Fig.1 The global coordinate system and co-rotationalcoordinate system of beam element

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

圖2 總體坐標(biāo)系與梁單元端截面隨動坐標(biāo)系Fig.2 The global coordinate system and co-rotationalcoordinate systems of the cross section of beam element

1.2 單元變形后隨動坐標(biāo)系的計算

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

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

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

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

根據(jù)圖3中各個矢量的關(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)單元端截面隨動坐標(biāo)系的計算

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

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

(2)

則各節(jié)點的轉(zhuǎn)動軸矢量方向為：

(3)

(3)單元隨動坐標(biāo)系的計算

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

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

(4)

(5)

(6)

1.3 單元隨動坐標(biāo)系中單元變形的計算

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

(1)彎曲變形

(7)

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

(3)伸長變形

在潘永仁的研究中，將單元的伸長變形定義為薄膜變形，通過定義單元變形曲線并計算曲線弧長變化來得到單元的薄膜變形。

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

已知t時刻單元節(jié)點位置為tXji，則可以計算得到t時刻單元長度為

(8)

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

Δl=t+Δtl-tl

(9)

Δl即為t+Δt時刻單元的伸長值。

1.4 節(jié)點位移計算

(1)單元節(jié)點力的計算

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

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

(10)

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

fe=keue

(11)

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

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

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

(3)節(jié)點位移計算

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

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

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

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

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

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

R(u)=F(u)

(12)

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

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

2.1 基本理論

在3維桿、梁單元非線性靜力計算的基礎(chǔ)上，本文開發(fā)了基于CR列式的非線性Newmark-β有限元動力時程積分方法，在算法中也考慮了結(jié)構(gòu)初始應(yīng)力的影響。其基本假定與普通的Newmark-β法相同：

(13a)

(13b)

(14)

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

(15)

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

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

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

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

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

K=a0M+a1C+KK

(16)

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

2.2 算法流程

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

圖6 非線性有限元動力時程計算流程圖Fig.6 Calculation flow graph of nonlinear time-history FEM

3 算例驗證

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

3.1 算例1 彎曲梁3維空間靜力計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 索-梁組合結(jié)構(gòu)面內(nèi)非線性振動ANSYS與本文程序計算結(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é)點力后結(jié)構(gòu)發(fā)生了振動，拉索在梁的帶動下發(fā)生了1∶1主共振，拉索端部(節(jié)點5)沿整體坐標(biāo)系Y方向0.01 m幅值的位移激勵，導(dǎo)致拉索中點(節(jié)點11)沿局部坐標(biāo)系Y1方向發(fā)生了0.06 m的振動。結(jié)構(gòu)振動體現(xiàn)了“拍振[17]”的非線性振動性質(zhì)，振動能量在拉索與整體結(jié)構(gòu)之間互相傳遞，呈現(xiàn)“此消彼長”的趨勢。符合非線性振動的理論預(yù)期。

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

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

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

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

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

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

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

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

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

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

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

3.3 算例3:拉索3維空間回旋運動計算

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

圖15 拉索受端點位移激勵示意圖Fig.15 The schematic figure of the cablewith displacement stimulation at the end point

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

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

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

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

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

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

造成上述現(xiàn)象的原因為：由于拉索的索力較大，在自重作用下，拉索面外1階自振頻率仍然接近面內(nèi)1階自振頻率，而拉索端點位移激勵沿圓周運動且頻率接近拉索的面內(nèi)與面外1階自振頻率，因此，在端點位移激勵下拉索發(fā)生了1階面內(nèi)-面外耦合的1∶1主共振，導(dǎo)致了較大幅度的3空間回旋運動，即：柔性索的“跳繩”運動現(xiàn)象。

4 結(jié) 語

(1)開發(fā)了3維空間桿梁單元基于CR列式的非線性Newmark-β動力時程有限元算法，闡述了程序計算原理，編制有限元計算程序，通過算例與ANSYS計算結(jié)果對比，驗證了程序的正確性與可靠性。

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

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

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