基于有限元-向量式有限元的斜拉橋非線性振動計算方法

王 濤， 胡宇鵬， 張興標， 劉德貴

(1.西南科技大學 土木工程與建筑學院，四川 綿陽 621000；2.中國工程物理研究院總體工程研究所，四川 綿陽 621999)

考慮幾何非線性的數值計算方法在實際工程的分析中有重要的應用。陳政清等[1]詳細研究了基于UL列式的非線性有限元(finite element,FE)法并將其運用到了橋梁工程的計算中。為了更簡便地進行非線性有限元算法編程開發，依據Wempner[2]提出的隨動坐標系格式方法(co-rotational formulation，CR列式)，鄧繼華等[3]開發了基于CR列式并考慮徐變作用的非線性桿、梁單元計算程序，應用到實際工程的靜力計算中。 潘永仁[4]使用CR列式幾何非線性算法開發了有限元計算程序，運用到大跨度懸索橋施工監控的靜力計算中。

王濤等[5]研究開發了基于CR列式的動力非線性有限元方法，動力時程求解中使用了在Newmark-β算法下的Newton-Raphson隱式動力平衡迭代。該方法求解穩定性好，可以使用較大的時間積分步長，可以計算柔性結構的大幅度非線性振動。非線性動力有限元法在動力時程計算的每一步中都需要計算總體剛度、質量、阻尼矩陣集成切線剛度矩陣。由于隱式求解要求矩陣不能奇異，所以通常不能直接計算結構斷裂、倒塌，破壞等無約束狀態。

Ting等[6]提出了向量式有限元(vector form finite element，VFFE)的概念，丁承先等[7]詳細論述了向量式有限元的特點。向量式有限元將結構視為在空間中通過單元聯結成為一群質點集合，單元無質量，只視為質點相互關系的依據，質點運動滿足牛頓第二定律。各個質點的內力、外力、加速度力依據單元聯結信息分別集成為一組向量，使用顯式的中心差分法得到結構動力計算結果。其優點在于，單元節點排序與計算規模無關。構件可以不限為桿件，可以處理大變形與大位移。在計算中可以改變結構單元的性質。由于不集成總體矩陣，可以在任意時刻切斷單元聯系，獲得結構的斷裂與倒塌行為。

關于向量式有限元，近幾年各個學者進一步發展了其在結構工程中的應用。陳沖[8]使用向量式有限元模擬計算了大型網架結構的倒塌與破壞。王震等[9-10]基于向量式有限元思想，開發了計算程序，使用薄膜單元與四面體實體單元研究了結構的非線性受力行為以及破壞與穿透狀態。曲激婷等[11]基于向量式有限元方法研究了黏滯阻尼減震結構抗豎向連續倒塌。Duan等[12]使用向量式有限元研究了列車-橋梁的耦合振動效應。李效民等[13]使用向量式有限元研究了淺海豎L型懸空管道三維動力響應。

離散元方法使用了更為一般的概念[14]，認為分析的結構對象是由一系列質點構成，質點之間的相互作用由接觸力確定，在巖土工程中應用較多，可用于求解土體顆粒之間的相互作用，模擬堆積、塌方等工程現象。為了在結構工程中應用離散元方法，齊念等[15-16]建立了桿、梁結構的離散元，認為質點分布于單元節點上，質點之間的接觸力以結構力學中桿、梁單元的計算公式來定義。在桿、梁結構的理論體系下，本文認為離散元與向量式有限元在算法本質上是一致的，同樣是針對質點的求解，不需要組集總體剛度矩陣，使用顯式數值積分計算結構動力問題，適合并行計算，葉繼紅等[17]開發了使用桿、梁離散元的GPU并行計算方法。

綜上所述，本文認為，使用桿、梁力學模型的向量式有限元、離散元，結構構件的內力與外力計算方法上是一致的，是同一個方法的不同概念。所以，在后文中本文只討論CR列式非線性有限元與向量式有限元的相同點與差別。

由于不組集總體剛度、質量矩陣，向量式有限元不能對結構進行模態計算。使用緩慢加載的顯式動力積分來近似處理靜力狀態，得到結構自質量下的靜力狀態也較為麻煩。

本文使用了新的編程開發思路，基于核心原理的一致性，將有限元、向量式有限元納入統一的程序框架。這樣，避免了如丁承先等研究中專門針對某個問題編程的重復工作。使用通用的節點、單元、實常數、邊界條件等的描述方法建立計算模型。使用非線性有限元方法計算結構的靜力狀態、動力特性，使用平衡迭代的有限元Newmark-β隱式方法計算結構長時間的非線性振動響應。通過直接設置調用向量式有限元求解模塊，即可模擬計算結構的沖擊、破壞、斷裂等狀態。本文方法彌補了有限元與向量式有限元各自的缺點，具有更好的適用性。

1 向量式有限元

如圖1所示，一根梁單元在外力作用下，在t1時刻，單元1、2節點發生了位移和轉動，經過Δt到達t2時刻變形狀態。依據向量式有限元虛擬逆向運動的思想，假定X2Y2坐標系平移、旋轉回到X1Y1坐標系的位置，兩個坐標系重合，則X2Y2中梁發生的變化為單元的純變形，這與CR列式非線性有限元方法通過扣除剛體轉動得到單元的實際變形狀態想法是完全一致的。單元節點位置即為向量式有限元的質點。

(a)

使用向量式有限元的求解規則如下。

(1) 梁的純變形。為了求解梁的純變形，令t1時刻單元做逆向的平移與轉動，使t1時刻單元的節點1與t0時刻單元重合，X2Y2與X1Y1方向重合。忽略單元的弦向變形，單元的伸長使用軸向變形來描述

Δl=l1-l0

(1)

單元逆向轉動的轉角為

Δφ=φ2-φ1

(2)

可以得到單元的實際轉角為

(3)

t1與t2時刻轉角的增量

(4)

(2) 單元內力計算。向量式有限元將運動歷程劃分為一系列途徑單元。在同一途徑單元內，內力計算以本時刻t1與下一個時刻t2構型為基礎。依據結構力學原理可以求得各個節點在局部坐標系下的內力增量。

(5)

式中,E、A、I為單元在t2時刻的材料切線模量、單元面積、為單元截面慣性矩。依據向量式有限元的思想，內力增量在結構變形的任意時刻t1都可以變化的。

考慮單元靜力平衡條件可以計算出節點內力為

ΣFx=0,f1x=-f2x

(6)

ΣFy=0,f1y=-f2y

使用坐標轉換可以得到全局坐標系下的單元的節點內力即

(7)

式中:Fi為總體坐標系下的節點內力向量；R為t2時刻總體坐標系到單元局部坐標系的坐標轉換矩陣；如圖1所示，傾角為φ2，所以

(8)

(3) 運動公式。 根據牛頓第二定律，點i的運動公式如下

(9)

(4) 求解方法。使用向量式有限元不需要組集總體剛度矩陣，通過結構模型單元連接信息集成總體內力向量F、外力向量P、質量向量m、位移向量x。向量式有限元一般使用中心差分法顯式求解，當考慮阻尼時，由中心差分公式得到

xn=Clh2m-1F+2C1xn-1-C2xn-2

(10)

式中:C1=1/(1+0.5ξh);C2=C2/(1-0.5ξh);ξ為阻尼比;h為積分時間步長。由式(1)～式(9)可知，向量式有限元直接建立了結構的外力向量、內力向量、加速度力向量的對應關系，使用顯式動力學積分，不求解矩陣方程，不存在有限元方法中矩陣奇異問題。向量式有限元求解規模與節點數量相關，與節點排序無關。

(5) 向量式有限元求解基本流程,如圖2所示。

圖2 向量式有限元計算流程

(6) 斷裂模式。在丁承先等、陳沖、王震等的研究中將斷裂模式設置為在模型節點處斷開，單元數量不變，但這會導致計算中臨時增加節點，對編程不利?？梢詤⒖紬U梁結構離散元模式將單元劃分更為密集，當單元力大于或定義某個時刻的時候，直接把單元力置為零，即可模擬斷裂情況。這時在單元處斷開不需要隨時間增加節點力向量，雖然增加了計算量，但能更符合實際情況。

(7) 碰撞檢測。結構斷裂后的構件可以自由運動，可能與其他構件或掉落地面發生碰撞。碰撞須滿足動量守恒與能量守恒定律。質點碰撞后對當前時間步質點的速度進行修正計算，在后續計算中若不發生碰撞，繼續使用式(10)計算，若檢測到碰撞則繼續修正速度。具體方法可參考陳沖的研究。

2 基于CR列式的非線性隱式有限元動力時程計算方法

依據參考潘永仁、王濤等的研究，本文作者開發了基于CR列式的非線性靜力、動力有限元算法。對于非線性有限元靜力計算，結構的總體的計算公式為

F(x)=P(x)

(11)

式中:F為結構內力向量;P為結構外力向量;它們都與節點位移x相關?？梢允褂脝卧獎偠染仃嚱M集結構總體切線剛度矩陣，使用Newton-Raphson方法進行平衡迭代得到結構的靜力平衡狀態。

對于每個二維梁單元，單元坐標系稱為隨動坐標，根據單元位置實時計算，單元的實際內力計算與向量式有限元思想相同，通過扣除剛體位移，得到單元實際位移向量。

xe={0 0θ1u120θ2}T

(12)

式中:θ為扣除剛體位移后，單元坐標系下實際轉角;u12為單元的伸長。計算方式與向量式有限元相同，依據有限元法規則，與式(5)的力學意義相同，使用單元剛度矩陣ke直接乘以單元位移向量

fe=kexe

(13)

得到單元的內力向量fe。與式(7)相同，依據當前時刻單元位置，使用單元隨動坐標系下的坐標轉換矩陣Te可得到總體坐標系下結構各個節點內力向量F與外力向量P。關于二維梁單元的剛度矩陣、坐標轉換矩陣等計算原理可以參考有限元方法文獻[18]。

對于動力求解問題，參考王濤等研究所述的詳細研究開發過程?；贑R列式考慮幾何非線性時，由t1時刻經過Δt結構t2時刻的振動方程為

(14)

式中：M為總體質量矩陣；C為總體阻尼矩陣；F為結構內力向量；P為結構外力向量。依據Newmark-β法設定t1、t2時刻位移、速度、加速度之間的關系為

(15)

(16)

式中，a0～a5為積分參數，將式(15)代入式(14)得到

(a0M+a1C)x(t2)+F[x(t2)]=P[x(t2)]+

(17)

在每一個積分時間步內求解非線性方程式(17)即可得到結構的非線性振動時程，式(17)左右端都與各個節點當前位移向量x(t)相關。式(17)左端記為結構等效內力向量F，右端記為結構等效外力向量P。在計算結構總體切線剛度矩陣后，即可使用Newton-Raphson法迭代求解。

根據以上論述，可知依據CR列式的靜、動力非線性有限元方法在結構的單元節點內力與外力以及坐標轉換的計算方法與向量式有限元完全一致。不同點在于，CR列式非線性有限元有限元法中，靜力使用內外力平衡迭代求解，動力是在Newmark-β算法框架下進行平衡迭代，為有限元隱式方法。而向量式有限元本質上是基于單元隨動坐標系考慮幾何非線性的顯式有限元動力時程積分。

3 有限元-向量式有限元統一算法框架

CR列式非線性Newmark-β隱式方法具有較好的收斂性與計算穩定性，可以使用較長的時間步長，可以較好地模擬結構低頻非線性振動響應。而向量式有限元不集成總體剛度矩陣，使用顯式動力學積分計算，通常使用較小的時間步長，可以獲得結構在短時間內受到沖擊、斷裂狀態下的動力響應。本文開發的有限元-向量式有限元統一程序框架,如圖3所示。

如圖3所示，向量式有限元在計算單元內力與外力時可以與CR列式的非線性有限元動力時程積分共用計算模塊。使用有限元-向量式有限元統一的算法框架，無需如丁承先等研究中那樣針對每個結構問題進行單獨地向量分析與建模編程。程序中只需建立有限元幾何模型，然后根據計算需要選擇進入向量式有限元計算流程即可。

在結構的動力分析中，通常需要首先得到結構在自質量下的靜力狀態。向量式有限元是不組集總體剛度、質量矩陣的顯式動力計算，計算結構靜力狀態時，需設定適當的耗散阻尼，使用緩慢加載來模擬靜力過程，過程相對較為繁復，也無法進行模態計算。所以，如圖3所示，本文使用該程序的分析流程為：首先，使用非線性有限元計算得到結構的靜力狀態；然后計算結構模態；最后再根據分析對象的實際情況選擇使用有限元或向量式有限元計算結構的非線性振動狀態。

向量式有限元為顯式積分方法，計算過程中無隱式積分那樣使用平衡迭代在每一步迭代中檢查收斂性，即：每一步均檢查式(17)中左端等效內力與右端等效外力是否達到平衡狀態。因此，算法的嚴謹性相對較為欠缺，時間步長的選取對計算結果影響較為明顯?？梢允褂秒[式方法來驗證向量式有限元求解時間步長的可靠性。

4 算例驗證

依據第1～第3章所述原理開發了有限元-向量式有限元的相關算例，為了與向量式有限元保持一致，有限元計算中使用集中質量矩陣。

4.1 彈簧擺自由振動斷裂后自由落體以及碰撞(算例1)

圖4 彈簧擺自由振動模型

圖3 有限元-向量式有限元統一算法框架流程圖

Fig.3 The flow chart of finite element & vector form finite element method unified algorithm framework

彈簧擺使用桿單元模擬，質量分布在節點上。單元彈性模量E=100 Pa，截面積A=1.0 m2，質量密度ρ=1.0 kg/m3，重力加速度G=9.8 m/s2，單元初始軸力為50 N，向量式有限元設置阻尼比ξ=0.2，計算時間步長h=0.000 1 s，求解50 000步，假設斷裂部分自由落體與地面接觸后碰撞豎向的速度全部反彈。

使用本文程序建立有限元模型，首先，使用非線性有限元靜力計算程序得到結構的初始構型，然后使用向量式有限元計算模塊得到結構的運動軌跡(動畫)如圖5所示。

圖5 彈簧擺自由振動的運動軌跡

Fig.5 The time-history of free vibration of the spring pendulum

從圖5可知，計算中彈簧擺首先使用有限元法得到了靜力狀態初始構型，釋放后發生了非線性自由振動，當單元4斷裂后，端部一個單元自由落體與地面發生了碰撞并反彈，圖中虛線為節點6運動的軌跡線。

得到節點6以及節點4的動力時程響應曲線,如圖6所示。

由圖5與圖6可知，向量式有限元較好地模擬了彈簧擺的運動過程。

(a) Y方向位移時程曲線

4.2 索梁組合結構非線性振動(算例2)

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

(2) 從實驗和數值角度研究了0Cr18Ni9不銹鋼材料實際的起裂斷裂韌度，針對I型裂紋緊湊拉伸測試，利用掃描電鏡斷面觀察和聲發射技術確定起裂載荷為23kN，根據不同方法所計算的J積分存在差異性。按照GB/T 21143標準計算的J積分的實際起裂斷裂韌度值JiGB為351.4kJ/m2，處于有限元計算的J積分范圍內。利用數字圖像相關方法測量的位移場計算得到的實際起裂韌度值JiDic(293.4kJ/m2)與有限元計算的J積分平均值JiFem(293.1kJ/m2)接近。

圖7 索-梁組合結構有限元-向量式有限元公用模型(m)

總體坐標系為XY，建立拉索模型局部坐標系為X1Y1。梁分為5個二維梁單元，拉索共分為20個二維直桿單元。不計泊松比，梁單元的物理參數為，彈性模量E=2.0×1011Pa，剪切模量G=1.0×1011Pa，材料質量密度ρ=7 800 kg/m3，單元截面積A=0.045 m2，抗彎慣性矩Iz=1.7×10-3m4。桿單元彈性模量、質量密度與梁單元相同，截面積A=3.0×10-4m2初始軸力H=5.0×104N。

向量式有限元算法無法計算結構自振頻率。所以，首先，使用非線性有限元計算結構的靜力狀態。然后，對結構進行模態計算。結構自振頻率與振型,如圖8所示。

(a) 第1階頻率:3.135 9 Hz

本算例對拉索的索力以及結構的總體剛度進行了設計，所以，模態計算結果表明，結構的1階、2階模態中可以看到拉索和整體結構同時發生振動的振型。依據Nayfeh等[19]的研究可以知道，在外激勵作用下，拉索會在梁的帶動下發生1階共振。

先進行靜力計算，然后在圖6的梁端點處施加Y方向外力-10 kN，然后釋放，使用本文開發的非線性有限元-向量式有限元計算得到拉索1/2點振動時程曲線如圖9所示。其中非線性有限元Newmark-β動力時程積分步長為0.001 s，向量式有限元非線性顯式動力學積分時間步長為0.000 1 s。

(a) 梁端點振動時程圖

由圖9可知，梁端受到集中力變形后，索-梁組合結構發生了自由振動，由于拉索局部頻率與整體結構頻率接近。振動能量在拉索和梁之間傳遞，梁與拉索的振幅呈現“此消彼長”的狀態。由頻譜圖可以看出，梁的結構的響應頻率3.15 Hz幅值較大8.60 Hz幅值較大，索的主要響應頻率為3.20 Hz。與計算得到的自振頻率接近，這說明，由于結構振幅較小，為弱非線性振動，接近線性振動響應。有限元與向量式有限元計算結果一致。

為了對比結構在較強非線性振動狀態下，兩種計算方法的計算結果。不設置阻尼，在圖7模型梁右端部作用沿y方向的豎向簡諧力P=P0sin(ωt)，其中外力項P0分別取-1.0 kN、-2.0 kN、-3.0 kN，激勵頻率ω取3.16 Hz。結果如圖10所示。

(a) 振動時程圖(P0=-1.0 kN)

由圖10可知，由于外激勵與整體結構以及拉索的頻率接近(3.136 Hz)，結構發生了共振。雖然未設置阻尼，但由于非線性振動效應導致的“振動硬化”現象，拉索的振幅不會一直增加，而是呈現了“拍振”的狀態，由圖10(b)可以看出，拉索發生了非線性振動，結構的響應頻率發生了改變，變得更加復雜，結構響應頻率的復雜度與振幅相關。關于非線性振動現象的更多理論解釋，可以參考Nayfeh等的經典著作。

觀察圖10(a)可以看出，分別使用有限元與向量式有限元計算，振動發展的初期，兩者計算結果時程曲線是重合的。但隨著振動的發展，差別相對變得明顯，依據Nayfeh等的研究，這是由于非線性振動對初始條件敏感決定的，微小的初始條件差別會隨著振動發展非線性增加。本文有限元、向量式有限元結構內力、外力計算方法相同，但由于動力算法的差別(隱式、顯式算法)，導致了振動時程圖差別隨時間發展逐漸明顯。

從各個頻譜圖可以看出，兩種計算結果的頻率響應趨勢是接近的，都呈現了非線性振動導致的響應頻率復雜化狀態。外力P0取值越大，頻譜圖中得到的結構頻率響應越復雜，振動非線性的程度越大。由時程圖可以看出，當P0分別為-1.0 kN、-2.0 kN、-3.0 kN時，時程曲線開始產生差別的時間分別約為第15秒、10秒、4秒。這說明非線性程度越大，非線性振動造成的差別放大效應越明顯，這符合非線性動力系統的特性。

綜上所述，本算例使用向量式有限元得到了結構非線性振動的結果，與CR列式非線性有限元對比，驗證了向量式有限元法的合理性與可靠性。

5 基于有限元-向量式有限元的車-橋耦合動力作用下斜拉橋非線性振動分析

纜索承重結構的斷索狀態分析是近期工程結構研究中的熱點內容。張超等[20]研究了多重四邊形環索-張弦穹頂局部斷索沖擊，對于拉索斷裂模擬，使用了單元退出工作的等效荷載瞬時卸載法。張羽等[21]研究了大跨混凝土斜拉橋施工過程中結構的斷索后的受力狀態，基于ABAQUS有限元軟件，使用靜力計算后瞬時拆除構件方法進行了研究。Zhou等[22]使用結構模型單元狀態迭代更新的等效荷載法，基于SAP2000有限元程序，開發了復雜的算法流程，分析了汽車作用下拉索斷裂時斜拉橋的非線性振動狀態。

以上研究中均采用了商業有限元程序進行研究，要得到合理的計算結果，技術上處理較為復雜。而使用本研究開發的有限元-向量式有限元統一算法框架，由于向量式有限元使用顯式動力時程計算的特性，可直接進行動態模擬，不需要進行復雜的幾何模型處理、受力邊界條件處理。只需要設定拉索斷裂的時間或受力條件，程序可以自動運行，達到條件即將單元內力置零來模擬拉索斷裂，可以方便地直接計算得到斜拉橋拉索斷裂時的非線性振動狀態，同時得到斷裂拉索在重力作用下的自由墜落狀態及其與橋梁結構可能發生的碰撞。

參考李永樂[23]所述列車-橋梁耦合動力作用計算方法以及秦順全[24]研究中的研究數據，以天興洲長江大橋為研究對象，建立二維簡化計算模型,如圖11所示。

圖11 天興洲長江大橋二維模型(m)

總體坐標系為XY，X為縱向，Y為豎向。斜拉橋為3索面，所以二維模型參數按照1/3簡化，各個拉索的編號如圖11所示。橋塔上靠近中跨、邊跨最長的3根拉索分為40段，其他拉索分為20段。 全橋共1 221個節點，1 433個單元。斜拉橋主梁端部豎向約束，主梁與輔助墩為豎向約束連接，主梁與橋塔橫梁位置為豎向約束連接，縱向使用彈性阻尼連接模擬縱向液壓阻尼器的約束。重力加速度取G=9.8 m/s2。阻尼比設為0.05%。計算列車使用二維4軸模型如圖12所示。

圖12 一節列車車廂模型

高速列車的剛度、阻尼、質量參數使用CHR2動車組，一列列車共8節車廂，與橋梁相同，也按照1/3簡化。軌道不平順使用德國譜低干擾譜模擬。使用車、橋系統分離迭代方法計算列車、橋梁相互作用力。

使用非線性有限元計算結構在自質量作用下的靜力構型，再計算得到全橋前2階自振頻率如圖13所示。

由圖13可知，第1階為全橋振動模態，第2階為拉索局部振動模態。全橋模型與拉索的局部頻率有明顯差別，根據王濤[25]的研究，在日常運營中，列車作用下斜拉橋不會在端部位移激勵下發生共振。

(a)

使用線性、非線有限元、向量式有限元分別計算一列列車過橋，計算結果如圖14所示。其中，列車速度為200 km/h，列車第1輪對進入橋梁開始計算，總共運行1 400 m。有限元計算使用空間步長0.25 m換算得到時間0.004 5 s，計算步數為5 600步。向量式有限元時間步長取0.000 1 s，列車運行距離換算為計算步數為252 000步。

(a) 主梁1/2點Y方向振動時程

由圖14可知，分別使用3種計算方法，得到的計算結果差別很小，其中線性、非線性有限元計算結果有微小差別，非線性有限元與向量式有限元計算結果一致，這是由于列車過橋時橋梁受力接近一個靜力過程，且非線性程度很小的原因。

依據前期各個研究，斜拉橋拉索受到火災或疲勞破壞通常更容易在端部斷裂。所以，這里在計算斜拉橋靜力狀態后，假設拉索32靠近主梁端部附近突然發生斷裂，使用向量式有限元計算，計算中不考慮拉索斷裂后的碰撞。在計算開始1 s后拉索發生斷裂，向量式有限元計算時間步長為0.000 1 s?？偣灿嬎?5 s得到計算結果，如圖15所示。

(a) 主梁1/2點X方向振動時程圖

由圖15(a)可知，主梁發生了明顯縱向(X方向)振動，振幅在0.01 m以下，由于跨中左側最長拉索發生斷裂，所以振動的平衡位置為正方向。主梁與橋塔之間存在豎向支撐同時設置了縱向彈性液壓阻尼器，所以縱向響應幅值較小，且隨時間衰減。主梁1/2點豎向(Y方向)振幅約為0.2 m。由圖15(c)可知，拉索32在重力作用下發生了大幅度非線性下墜振動(不計碰撞)。由圖15(d)可知拉索第6秒的下墜狀態。

為了研究拉索斷裂對高速列車行駛安全性的影響。假定拉索斷裂時列車仍然沿軌道在橋面上行駛，當列車作用以200 km/h速度通過斜拉橋，使主梁1/2點位移達到最大時，拉索32靠主梁端部突然發生斷裂，得到計算結果如圖16所示。

由圖16(a)可知，拉索32斷裂時，斜拉橋主梁1/2點發生了明顯的較大幅度豎向振動，由于列車作用在橋梁上，振幅較圖15(b)更大。

由于這里假設列車始終在橋梁的軌道上行駛，觀察圖16(b)可知，當拉索發生斷裂時，列車的豎向振動幅值突然最大增加約0.2 m。觀察圖16(c)，列車的豎向加速度由約0.4 m/s2突然增大至接近4.0 m/s2。依據李永樂提供的研究數據，列車過橋的豎向加速度需小于3.63 m/s2。顯然，列車車廂豎向的位移、加速度值的突變會使列車乘客產生較大恐慌，拉索32的突然斷裂會對列車的行駛安全性造成威脅。

(a) 主梁1/2點Y方向振動時程圖

得到列車過橋拉索斷裂時各個拉索的索力變化結果,如圖17所示。

(a) 拉索31應力時程圖

由圖17可知，當高速列車過橋運行至跨中時若拉索32發生斷裂，拉索31的應力幅值相對較大，而右半跨的跨中拉索33以及左半跨跨中靠橋塔拉索17的應力幅值相對較小，但拉索承受的最大應力的均遠小于拉索設計極限強度1 860 MPa。

這說明，對于斜拉橋這個較為復雜的多自由度結構體系，少量拉索斷裂后，其他拉索能夠有效分散承擔結構受到的外部荷載與結構自質量，離斷裂拉索越遠的拉索受到的影響相對越小，所以，本文認為，列車過橋拉索32斷裂時，橋梁結構是相對較為安全的。

6 結 論

(1) 基于隨動坐標系(CR列式)的動力非線性Newmark-β有限元法為使用平衡迭代的隱式動力時程積分方法。而向量式有限元法本質上是CR列式幾何非線性有限元顯式動力時程積分方法。

(2) 有限元通常無法直接計算結構受到破壞、斷裂、無約束等非線性動力狀態，向量式有限元可以模擬結構的斷裂、破壞、碰撞等非線性動力響應，但無法進行直接靜力計算、結構模態計算?；诤诵脑淼囊恢滦?，本文開發了有限元-向量式有限元統一算法框架，彌補了這兩種方法各自的缺點。

(3) 對于大跨度斜拉橋，使用本文的有限元-向量式有限元結合的方法能更好地對結構進行全面的計算分析。對于列車-橋梁耦合振動，使用非線性有限元與向量式有限元的計算結果是一致的。本文研究表明：當高速列車運行通過斜拉橋時，若跨中最長拉索突然發生斷裂，斜拉橋主梁會發生較大幅度振動，對橋梁結構的安全性影響相對不大，但會對列車的行駛安全性造成威脅。

(4) 向量式有限元可以直接計算斜拉橋在極端情況下的非線性振動。計劃在后續研究中進一步開發使用桿、梁單元的三維空間向量式有限元算法程序，針對實際斜拉橋增加更多極端作用下的計算工況，為實際斜拉橋在動力破壞狀態下的安全性評估提供研究依據。