基于MATLAB的公路橋梁車橋耦合數值計算方法

易晉生，顧安邦，王小松

(重慶交通大學土木建筑學院，重慶400074)

車橋耦合振動的研究最早可追溯到19世紀初期。1847年當列車通過英國Chester鐵路橋時由于振動產生撓度過大而導致橋垮塌，此后100多年來，車橋動力相互作用問題一直是各國橋梁設計工作者研究的重點課題［1］。當車輛通過橋梁時，在橋梁自身慣性力和車輛自振慣性力等動載作用下，橋梁變形和受力會大于相同荷載靜力作用。尤其當車輛的激振頻率和橋梁自振頻率相接近時，會發生共振現象，使得橋梁產生過大撓度甚至破壞。

近年來，隨著橋梁計算理論和施工技術的不斷完善，以及高強輕質材料在橋梁工程中的應用，促使橋梁結構向著跨度越來越大、質量越來越輕、剛度越來越小的體系發展。隨著人們對交通運輸的需要日益增長，不僅車輛的數量迅速的增長，而且車輛的行駛速汽車與眾多的服役期滿或損傷的橋梁承載能力不足之間的矛盾也日益突出［2］。此外，振動過大將嚴重影響到行車的舒適性，更為嚴重的還將影響到行車安全。有關研究表明，現行規范關于橋梁結構沖擊系數的取值是偏不安全的［3］。這兩方面因素使得橋梁承受的汽車活載以恒載的比例越來越大，車橋耦合振動效應越趨突出，已成橋梁設計計算重要因素之一。

筆者根據所推導出的車橋動力方程組，利用數學分析軟件MATLAB，結合兩種計算精度較高而實用的數值計算方法編制了車橋耦合振動計算程序。針對用NEWMARK法求解分離的車輛和橋梁運動方程組時，在每個時間步長內迭代計算直至橋梁響應平穩以提高計算精度。為求解復雜的車橋耦合振動問題提供了簡便、實用而且計算精度較高的數值計算方法。

1 車橋模型及振動方程的推導

1.1 車輛模型及車輛運動平衡方程

人們在對車橋相互作用問題的研究過程中建立了各種車橋分析模型，其中，梁的模型都是彈性連續體，而對簡化車輛荷載就具有從簡單到復雜的各種形式，從最早的不變的移動荷載、考慮慣性力的移動質量，到后面的由彈簧和阻尼器所連接的簧上質量模型以及豎向密貼和豎向非密貼現代車輛模型［4－5］。筆者采用兩自由度汽車模型(圖 1)，將車輛視為多剛體組成的動力系統，質量集中在車體上，車輪和車體之間的支撐由彈簧和阻尼器來模擬。此車輛共兩個自由度，它們分別是車體的沉浮、點頭。

圖1 兩自由度車輛模型Fig.1 Vehicle model with 2-DOF

假設車輛坐標系以圖1中所示為正方向，車體、車架和輪對均視為剛體，不考慮各部件的變形;彈簧和阻尼均為線性;車輛各部件在各自平衡位置附近做小位移的振動。利用達朗貝爾原理對圖1所示的兩自由度車輛的振動方程進行推導如下。

車體豎向運動平衡方程為:

車體點頭運動平衡方程為:

對式(1)、式(2)進行整理并寫成如下矩陣方程的形式:

式中:Mc為車輛的質量矩陣;Cc為車輛的阻尼矩陣;Kc為車輛的剛度矩陣;Fc為車輛運動方程的荷載項。具體如下:

1.2 橋梁運動平衡方程

橋梁是多自由度體系，其運動平衡方程為:

式中:Mb為橋梁的質量矩陣;Cb為橋梁的阻尼矩陣;Kb為橋梁的剛度矩陣;Fb為橋梁運動方程的荷載項。

分別取前、后輪對作為研究對象，可以得到Fb的表達形式，限于篇幅就不一一展開，這里直接給出車輛對橋梁作用的荷載Fb:

2 車橋耦合數值方法及計算原理

車橋耦合振動的分析方法可以分為時域法和頻域法，時域法又分為以下兩種方法［6-7］:

1)將車輛模型與橋梁模型的所有自由度通過輪軌關系耦合在一起，消除不獨立自由度，建立統一的系統運動方程組，進行同步求解。

2)將車橋系統以輪軌接觸為界，分為車橋兩個子系統，分別建立車橋的運動方程，兩者之間通過輪軌接觸處的位移協調條件與輪軌相互作用力的平衡關系，采用迭代法求解系統響應。

從計算車橋耦合振動的數值分析方法中作者選用了常用的精度較高而且計算速度較快的龍格-庫塔法和NEWMARK法［8］，并結合它們各自的計算特點編制了MATLAB計算程序，并在后面通過算例進行結果比較。

2.1 龍格-庫塔法

龍格-庫塔法［9］是單步法中的一種，它精確度較高，無需計算高階導數，就可以達到與泰勒級數方法同樣的精度。其基本思路是構造在某些點處的值的線性組合公式，使其按泰勒展開后與初值問題的解的泰勒展開式比較，以相同項系數相等以確定其中的參數。

筆者主要對4階龍格-庫塔法的求解思路進行說明。并選用4階龍格-庫塔法進行計算，采用MATLAB編制了相應的計算程序，對4自由度橋梁與橋梁耦合進行了計算分析。4階龍格-庫塔法的計算公式如下:

2.2 NEWMARK 法

2.2.1 NEWMARK 法計算原理

NEWMARK法［10-11］是結構動力分析中最常用的計算方法之一，1959年N.M.Newmark基于公式(11)、公式(12)發展了一類時間步進法:

參數β和γ定義了時間步內加速度的變化，并決定方法的穩定性與精度特征。對于γ=1/2和1/6≤β≤1/4的典型選擇，從包括精度的所有觀點來看都是令人滿意的。式(1)、式(2)與時間步結束時的平衡方程結合，提供了從i時刻已知的ui，˙ui，¨ui計算i+1時刻的ui+1，˙ui+1，¨ui+1的基礎。執行這些計算需要迭代，因為未知的¨ui+1出現在式子的右側。然而，對于線性體系，修正NEWMARK法的原始公式可以允許使用式(11)、式(12)求解時不迭代。單自由度問題采用NEWMARK法求解的具體步驟如下:

1)初始計算

選擇Δt

2)對每個時間步i進行計算

3)對下一個時間步進行循環。而i由i+1取代，對下一個時間步重復2)。

2.2.2 采用NEWMARK法計算車橋分離系統策略

將車橋系統以輪軌接觸為界分為車橋兩個子系統，分別建立車橋的運動方程，采用迭代法求解系統響應，同時兩個系統之間在輪軌接觸處滿足位移協調條件與輪軌相互作用力的平衡關系。基于MATLAB平臺，采用NEWMARK法進行車橋耦合數值計算分析的具體求解思路與過程如下:

1)在t時刻，提取當前橋梁的動力響應作為初始迭代值，并以此為基礎通過插值計算車輛輪對處響應的位移、速度和加速度;

2)根據前面求得的車輛輪對處響應的位移、速度求車輛受到的力的作用，并利用NEWMARK法計算t+dt時刻車輛的動力響應;

3)根據前面利用耦合條件計算的輪對響應的位移、速度和加速度和t+dt時刻車輛的動力響應，計算得到車輛對橋梁的作用力;

4)利用t時刻的橋梁響應采用NEWMARK法計算t+dt時刻的橋梁動力響應;

5)重復步驟1)～步驟4)，直至兩次相鄰計算的橋梁響應收斂，再進行下一個時間步計算。

步驟5)中先根據前后相鄰兩次計算得到的同一時刻的橋梁響應的差值向量，然后計算其向量的范數作為收斂性的判斷的依據。步驟5)的收斂性判斷是必要而且關鍵的，這將在后面的計算中得到證明。

3 算例分析

橋梁參數為:跨徑L=30 m，抗彎慣性矩 I=8.65 m4，彈性模量 E=2.943 ×1010N/m，線質量m=3.64 ×104kg/m，泊松比 ν=0.2。

車輛參數為:Mc=5.4 ×105kg，Jc=1.38 ×107N·m2，K1=4.135 ×107kg/m，C1=0 N·s/m，K2=4.135 ×107N/m，C2=0 N·s/m，S1=8.75 m，S2=8.75 m。

計算參數為:橋梁運行速度v=100 km/h，模態阻尼比ζ=0，時間積分步長為0.01 s。

分別采用基于統一方程求解的龍格-庫塔法和分離迭代法求解的NEWMARK法對上述算例求解。圖2～圖4分別為兩種方法下簡支梁跨中位置的位移響應、速度響應、加速度響應比較。為便于對比，同時給出了采用NEWMARK法而未迭代計算至橋梁振動平衡狀態下的結果。

圖2 橋梁跨中的位移響應比較Fig.2 Comparison of response of vertical displacement at Mid-span

圖3 橋梁跨中的速度響應比較Fig.3 Comparison of response of vertical velocity at Mid-span

圖4 橋梁跨中的加速度響應比較Fig.4 Comparison of response of vertical acceleration at Mid-span

由圖2～圖4可知:

1)由龍格-庫塔法和采NEWMARK法迭代計算車橋分離系統都得到了較精確的結果，兩者之間的符合度較高。

2)采用NEWMARK法迭代計算車橋分離系統時，在每個時間步內迭代收斂速度較快，并沒有因為平衡迭代損失太多時間，并且結果精度較高。

3)采用NEWMARK法計算車橋分離系統時，所得到的跨中撓度位移值與龍格-庫塔法的結果最大偏離僅為0.6%;而未進行迭代計算至橋梁振動平衡狀態時，其與龍格-庫塔法的最大偏離為5.1%。可見，采用NEWMARK法進行迭代計算的計算精度比未進行迭代計算至橋梁振動平衡狀態的算法要高。因此，在每個時間步內，迭代計算求解分離的車輛和橋梁運動方程組直至橋梁振動平衡狀態是必要且有意義的。

4 結語

應用達朗貝爾原理推導了兩自由度車輛模型和橋梁的運動平衡方程，選用了龍格-庫塔法和NEWMARK法兩種常用的、精確度較高的數值計算方法求解車橋耦合振動。在采用NEWMARK法迭代求解車橋分離雙系統方程組時，選擇合理的收斂條件作為迭代計算是否終止的依據。通過對比可知，兩種數值計算方法都能獲得較高的精度。

同時，比較了采用NEWMARK法分析車橋分離系統時，在每個時間步內，迭代計算車輛-橋梁運動方程組直至橋梁振動平衡狀態，和未迭代計算車輛-橋梁運動方程組至橋梁振動平衡狀態兩者的計算結果。研究表明，在每個時間步內迭代計算車輛和橋梁運動方程組直至橋梁振動平衡狀態是有必要，其計算精度比未進行迭代計算至橋梁振動平衡狀態的算法要高。

［1］ 夏禾.車輛與結構動力相互作用［M］.北京:科學出版社，2002:1-5.

［2］ 卜建清，婁國充，羅韶湘.汽車對橋梁沖擊作用分析［J］.振動與沖擊，2007，26(1):52-55.

Bu Jianqing，Lou Guochong，Luo Shaoxiang.Analysis on impact effects of bridge due to moving vehicle ［J］.Journal of Vibration and Shock，2007，26(1):52-55.

［3］ 丁南宏，林麗霞，孫迎秋.公路連拱橋在單車荷載下振動的理論和實驗研究［J］.蘭州交通大學學報，2005，24(3):28-32.

Ding Nanhong，Lin Lixia，Sun Yingqiu.Theoretical and experimental study of highway arcade bridge’s vibration under single vehicle load［J］.Journal of Lanzhou Jiaotong University:Natural Science，2005，24(3):28-32.

［4］ Green M F，Cebon D.Dynamic response of highway bridges to heavy vehicles loads:theory and experimental vibration［J］.Journal of Sound and Vibration，1994，170(1):51-78.

［5］ 曹雪琴，劉必勝，吳鵬賢.橋梁結構動力分析［M］.北京:中國鐵道出版社，1987:60-90.

［6］ 韓萬水.風-汽車-橋梁系統空間耦合振動研究［D］.上海:同濟大學，2006.

［7］ 李小珍，強士中.列車-橋梁耦合振動研究的現狀與發展趨勢［J］.鐵道學報，2002，24(5):113-120.Li Xiaozhen，Qiang Shizhong.State of the art review and trend of studies on vehicle-bridge interaction［J］.Journal of the China Railway Society，2002，24(5):113-120.

［8］ 沈火明，肖新標.求解車橋耦合振動問題的一種數值方法［J］.西南交通大學學報，2003，38(6):658-662.Shen Huoming，Xiao Xinbiao.Numerical method for vehicle-bridge coupling vibrations［J］.Journal of Southwest Jiaotong University，2003，38(6):658-662.

［9］ Steven C C，Raymond P C.工程數值方法［M］.6版.于艷華，傅效群，趙紅宇，等，譯.北京:清華大學出版社，2010:728-737.

［10］Anil K C.結構動力學理論及其在地震工程中的應用［M］.2版.謝禮立，呂大剛，譯.北京:高等教育出版社，2005:124-127.

［11］徐榮橋.結構分析的有限元法與MATLAB程序設計［M］.北京:人民交通出版社，2006:238-246.