基于流動坐標系的3維空間動力非線性有限元方法

王 濤， 劉德貴， 胡安杰

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

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

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

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

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

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

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

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

2 基于CR列式桿梁結構幾何非線性計算原理

2.1 三種坐標系統

根據潘永仁等人的論述，本文認為，基于CR列式非線性有限元算法的核心問題在于建立與有限元模型單元相關的3種坐標系統(即：結構建模總體坐標系、單元隨動坐標系、單元端截面隨動坐標系)以及處理它們之間的關系。由于直桿單元屬于梁單元的退化情況，所以這里以3維梁單元為例來闡述桿、梁單元中3種坐標系之間的關系。

(1)結構總體坐標系

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

(2)單元隨動坐標系

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

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

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

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

1.2 單元變形后隨動坐標系的計算

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

(1)歐拉有限轉動公式

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

圖3 有限轉動矢量示意圖Fig.3 The schematic figure of finite rotation vector

根據圖3中各個矢量的關系可以求出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)單元端截面隨動坐標系的計算

已知在t至t+Δt時刻，在結構坐標系中單元節點發生了轉動位移矢量Δθji，平動位移矢量ΔUji，j=1,2;i=1, 2, 3。

根據Δθji由數學定義可知轉動角度值為：

(2)

則各節點的轉動軸矢量方向為：

(3)

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

設t時刻在總體坐標系下，單元節點坐標為tXji，則可以得到t+Δt時刻節點坐標為：

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

(4)

(5)

(6)

1.3 單元隨動坐標系中單元變形的計算

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

(1)彎曲變形

(7)

(2)扭轉變形

(3)伸長變形

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

本文認為，對于一般的幾何非線性計算，根據小應變假設，從直觀上來看可以由節點的位置直接計算得到單元的伸長變形值。

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

(8)

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

Δl=t+Δtl-tl

(9)

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

1.4 節點位移計算

(1)單元節點力的計算

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

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

(10)

單元隨動坐標系中單元實際內力為：

fe=keue

(11)

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

(2)單元坐標轉換矩陣與應力剛度矩陣

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

(3)節點位移計算

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

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

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

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

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

如果結構處于最終的靜力平衡狀態，有如下關系：

R(u)=F(u)

(12)

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

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

2.1 基本理論

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

(13a)

(13b)

(14)

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

(15)

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

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

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

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

使用Newton-Raphson迭代法求解式(15)，當迭代未收斂時，式子左右兩端是不相等的，計算兩端的差值就可以得到迭代的不平衡力。在迭代過程中使用的切線剛度矩陣為Newmark-β法的等效總體剛度矩陣：

K=a0M+a1C+KK

(16)

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

2.2 算法流程

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

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

3 算例驗證

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

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

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

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

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

圖8中給出了P=300,P=600時，結構在靜力作用下發生大幅度變形后的構型。

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

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

表1 計算結果對比

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

3.2 算例2：索-梁組合結構3維空間非線性振動計算

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

圖9 索-梁組合結構有限元模型Fig.9 Cable-Beam structure FEM modal

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

表2 結構單元物理參數

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

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

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

圖10 結構前6階模態Fig.10 The first six modes of the structure

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

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

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

圖11 索-梁組合結構面內非線性振動ANSYS與本文程序計算結果Fig.11 The results of the Cable-Beam structure in-plane nonlinear vibration which calculated by ANSYS and the program of this thesis

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

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

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

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

這里如果使用線性動力時程計算，由于結構剛度不變，在結構發生共振時，理論上振幅會持續增加，直致小阻尼限制的上限，會遠大于非線性振動幅值上限，也不會發生非線性振動特有的振動硬化現象。所以，柔性結構的大幅振動必須考慮幾何非線性效應。

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

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

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

圖14 索-梁組合結構非線性振動時程狀態Fig.14 The time history states of nonlinearvibration of the cable-beam structure

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

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

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

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

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

圖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維直桿單元。

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

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

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

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

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

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

4 結 語

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

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

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

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