基于虛擬激勵法列車脫軌系數概率統計研究

周子驥，張 楠，孫琪凱

(北京交通大學土木建筑工程學院，北京 100044)

在車-橋耦合系統動力響應研究中，軌道不平順是引起耦合系統振動的重要激勵源。由于軌道不平順隨機函數是平穩Gauss 隨機過程，因此引起的車橋耦合振動也是典型的隨機過程[1]。確定性分析只針對給定參數得到的單一結果，難以體現車橋系統安全性和平穩性指標的概率分布規律，必然存在一定的風險，由此開啟了基于車橋隨機振動的橋梁安全性研究。

傳統車-橋隨機振動研究常采用的Monte Carlo法，通過諧波疊加等方法將軌道不平順功率譜生成多組不平順時域樣本，根據不同的時域樣本計算車-橋系統動力響應，再采用Monte Carlo 法統計計算結果。此方法是基于大樣本容量的隨機分析方法[2]。Chatterjee 等[3]根據橋面隨機不平順激勵樣本，采用基于Monte Carlo 法計算橋梁沖擊系數的均值和偏差。雖然Monte Carlo 法原理簡單，在車-橋隨機振動研究中被廣泛使用，但是由于其必須取足夠數量的樣本才能保證統計結果的可靠性，效率極低[4]。在此基礎上，學者們又進一步運用重要抽樣法、一次二階矩等方法研究橋梁結構可靠度[5?6]。同時，部分學者通過改進迭代計算方法來提高計算的精確性與時效性[7]。隨著有限元仿真的大力發展，學者們利用有限元模擬建立車-軌系統研究了該系統的隨機振動特性[8]。為進一步降低隨機振動分析的計算量，同時提高隨機振動分析的準確性，林家浩等[9]提出了隨機振動的虛擬激勵法。

基于虛擬激勵法，余志武等[2,10]、李小珍等[11?12]、趙巖等[13]、何旭輝等[14]等將其用于車-橋隨機系統，研究了列車荷載作用下車-橋耦合系統的隨機振動特性，并采用3 σ法則對橋梁響應上下限值進行估計，研究了橋梁的隨機動力特性。王立東等[15]提出了一種考慮軌道隨機不平順的列車地面振動預測方法，將土體簡化為2.5D 有限元匹配模型，研究了列車隨機振動對土體的影響。在此基礎上，部分研究亦在車橋系統分析中的同時考慮了其他類型的隨機荷載，分析了隨機荷載對橋梁動力響應極值的影響[16?17]。但是大多數研究主要針對橋梁結構安全性，對于列車運行安全指標的研究很少。在結構安全方面的研究，常用方法是假設概率分布函數得到安全界限值，然后開展討論[18]。對于車-橋系統而言，列車安全事故發生概率要高于因列車運行造成橋梁結構破壞概率，而且脫軌系數不服從正態分布，需采用新的概率分布函數進行統計研究。因此對于列車行車安全指標的研究有較大的意義。

基于上述研究，本文采用豎向密貼橫向線性蠕滑輪軌關系建立車-橋動力學模型，采用虛擬激勵法求得列車輪軌力功率譜；運用諧波疊加法，通過每一時刻三角級數分解求非平穩輪軌力功率譜的時域樣本；再由極值理論對列車脫軌系數進行安全統計；最后選取置信度99.73%的脫軌系數極值與多次時域樣本最大值進行比較，驗證使用極值理論研究脫軌系數的正確性。

1 車-橋耦合系統計算模型與輪軌關系

將車輛和橋梁看成兩個子系統，分別建立車輛和橋梁系統的模型。

根據拉格朗日方程可建立車輛系統動力學方程，列車運動方程為：

式中：Mv、Cv和Kv分別為車輛系統的質量矩陣、阻尼矩陣和剛度矩陣；Fv為車輛系統受到的車橋系統相互作用力；X、X˙ 、X¨分別是位移向量、速度向量和加速度向量。車輛單元位移向量的自由度排列順序為：

式中：下標c 表示車體；t1、t2 依次表示前、后轉向架；w1、w2 依次表示與前轉向架相連的兩輪對；w3、w4 依次表示為與后轉向架相連的兩輪對。

在車橋耦合系統相互作用研究中，為降低計算工作量，常采用振型疊加法建立橋梁動力方程。橋梁系統模型采用空間梁單元形式，每個單元節點有6 個自由度，分別為沿著x軸、y軸、z軸的平移和繞著x軸、y軸、z軸的轉動。本文不考慮軌道結構，假設軌道結構和橋梁沒有相對位移，輪軌力通過軌道直接作用于橋梁上。基于此可以得到橋梁運動方程為：

式中：Xb是橋梁廣義坐標向量；w是頻率向量；Φ為振型矩陣。

本文輪軌相互作用模型采用橫向線性蠕滑，垂向密貼的輪軌關系假定。輪軌垂向假定即假定車輛子系統坐標z方向，輪軌間無相對位移，此方向的相互作用力可以由輪軌間該方向的相對運動狀態確定；橫向線性蠕滑理論即滿足虛擬激勵法對線性系統的要求，又可考慮輪軌間橫向運動和相應作用力之間的線性關系。本課題組研究發現，對于我國常見的車輪半徑，線性蠕滑參數與非線性蠕滑參數誤差在0.13%～4.43%。因此采用橫向線性蠕滑理論可以滿足工程精度要求[19]。

2 基于虛擬激勵法的脫軌系數統計研究

2.1 車-橋耦合時變系統隨機振動分析

由第1 節可知，列車子系統和橋梁子系統通過輪軌接觸關系組成整體耦合時變系統，其耦合時變系統動力方程：

式中：角標v 和b 分別為車輛子系統和橋梁子系統；M、K和C分別為質量矩陣、剛度矩陣和阻尼矩陣；F為力向量。

列車橋梁耦合系統所受外荷載F(t)包括車輛自重引起的確定性動荷載Fg和軌道不平順引起的隨機性荷載Fr(t)。即：

盡管車橋子系統都是線性系統，但是由于列車在橋上的位置是隨時間不斷變化的，因此車橋系統具有非平穩特性。可將隨機激勵F?r(t)表示為：

式中：G(t)=diag[g(t?t1);g(t?t2);···;g(t?tn)]；x(t)=[x(t?t1),x(t?t2),···,x(t?tn)] ；Γ(t)為作用力指示向量，其作用是按輪距成反比的加權因子將輪對作用力分配到相鄰的兩個橋梁節點上；x(t)為平穩隨機過程；G(t)為調制函數矩陣；n為輪對總數。

假定u?i(t) 與u?j(t)是系統任意時刻響應向量，則系統隨機響應功率譜矩陣Suiuj(w,t)：

式中：hi(ti?τi,τi)為第i時刻脈沖響應函數。當i=j時，Suiuj(w,t)表示i時刻響應的自功率譜；當i≠j時，Suiuj(w,t)表示i、j時刻響應的互功率譜。

最后，根據隨機響應功率譜密度矩陣，通過下式求得系統隨機響應的均方根：

式中：Kz1、Cz1、Mw分別為車輛一系懸掛剛度、阻尼及輪對質量；Cc為輪軌間橫向蠕滑而產生的附加阻尼；Sa、Sv、Sc分別代表軌道方向、高低、水平不平順功率譜；b1為左右軌距中心距之半。橫向力與豎向力標準差：

式中： σz、 σy分別代表豎向力與橫向力的標準差；Szz、Syy分別表示由軌道高低不平順引起的豎向力功率譜和軌道方向不平順引起的橫向力功率譜；Suz、Suy表示由軌道水平不平順分別引起的豎向力功率譜和橫向力功率譜；dw頻率間隔帶寬。

2.2 非平穩功率譜的諧波疊加法

采用諧波疊加法模擬隨機過程是常用手段。土木工程領域常采用此方法模擬軌道、路面不平順及隨機風場。此方法要求功率譜函數是平穩過程[20]。對于輪軌力而言，由于橋梁位移變化具有時變特性，輪軌力功率譜也具有時變性，因此輪軌力功率譜函數，不能直接使用。

對于非平穩功率譜函數，引入調制函數，采用三角級數疊加法在每一時刻進行諧波疊加求其時域序列x(t)：

式中：x(t)為輪軌力時程序列；Sf(w)為輪軌力功率譜；g(t)調制函數；dw頻率間隔帶寬；wk為第k個對應頻率值； ?k為響應第k個頻率的相位，一般按照0～ 2π間均勻分布取值； σF為非平穩激勵下的標準差； σFs為平穩激勵下的標準差。

2.3 極值理論

由概率論，服從極值分布需同時滿足三個條件：觀測對象是隨機變量；這個隨機變量的底分布保持不變，或者如果有任何變化，可通過數據變換減少其影響；觀測到的極值是獨立的[21]。對于列車脫軌系數分布而言，滿足上述條件1 和條件3。對于高速鐵路橋梁而言，橋梁撓跨比值受到嚴格控制，因此橫、豎向輪軌力受軌道不平順的影響十分顯著，其分布高度近似服從高斯分布。脫軌系數是橫向力與動輪重的比值，其安全指標為[?0.8，0.8]，對于動輪重和橫向力，可以通過線性變換使其轉化為標準正態分布，因此二者比值服從柯西分布。對于脫軌系數而言，滿足極值分布的第二條件。由此可以得出，脫軌系數滿足極值分布條件。

式中：u是位置參數； σ是尺度參數； γ為歐拉常數，取0.5772；Dx為極值分布的方差；Ex為極值分布的均值。

極值分布函數求解流程圖，見圖1。

圖1 極值分布函數求解流程Fig. 1 Solution flow of extreme value distribution function

3 算例驗證及計算結果分析

本文以秦沈客運專線現場實驗為工程背景，通過與實驗結果比較來驗證本文模型與程序的準確性。因此，本文以秦沈客運專線狗河特大橋單跨24.6 m 標準簡支梁為例；計算列車采用先鋒號列車動車組，列車采用6 車編組，即動+拖+動+動+拖+動。運行時速270 km/h；線路橫向偏心2.5 m，垂向偏心1 m，其輪軌力與橋梁作用位置關系圖如圖2 所示；輪軌間激擾采用由德國高速鐵路低干擾譜變換出的軌道不平順時域樣本。車輛上橋前，先以與橋上相同的線路條件行駛150 m，待車體振動趨于穩定后進入橋跨結構，出橋后車輛再行駛450 m。橋梁與車輛參數見文獻[19]。

圖2 作用在橋梁上的力向量示意圖 /cmFig. 2 Schematic diagram of force vector acting on bridge

由圖3、圖4 為本文計算結果和文獻計算結果，二者比較可以看出，橋梁跨中豎向動位移幅值和變化趨勢基本吻合。跨中豎向加速度變化趨勢相近，幅值與文獻計算結果有差異，主要原因是由于文獻采用的是實測軌道不平順，而本文采取德國高速鐵路低干擾譜，軌道不平順對橋梁加速度有較大影響，因此導致本文計算結果和文獻結果稍有差異。通過上述計算結果與文獻結果對比可驗證模型與程序的正確性。

圖3 跨中豎向動撓度Fig. 3 Vertical dynamic deflection in mid-span

圖4 跨中豎向加速度Fig. 4 Vertical acceleration in mid-span

基于上述車-橋模型研究橋梁動力響應的隨機特征。其中軌道不平順為德國低干擾軌道不平順譜，波長范圍為1 m～80 m，相應的頻率計算范圍為0.9375 Hz～75 Hz，頻域積分步長為0.9375 Hz，考慮輪對慣性力的影響。

3.1 橋梁變形對輪軌力平穩特性的影響

由輪軌密貼與橫向線性蠕滑理論可知，豎向輪軌力含w2項，而橫向輪軌力僅含w項，因此，高頻對豎向輪軌力有較大影響。考慮軌道不平順二階導數項，結果如圖5 所示。考慮軌道二階導數項時，虛擬激勵法求得的豎向力均方根值與100 次時域循環求得均方根值較為吻合。因此，采用虛擬激勵法研究考慮輪對慣性力時的車橋系統隨機振動，密貼模型能滿足要求。

圖5 首輪對豎向輪軌力均方根Fig. 5 Root mean square of vertical wheel rail force of first wheel pair

由圖6 可知：橫向輪軌力功率譜首次峰值出現在20 Hz 附近，頻率越高，譜值越小；豎向輪軌力功率譜主要集中在中高頻階段，首次峰值出現在8 Hz 附近，當頻率大于28 Hz，輪軌力譜值趨于穩定。對于密貼模型，輪軌力受高頻影響較大。首輪對經過橋梁時間段為0.4 s～0.728 s，此段時間，橫、豎向輪軌力功率譜波動不明顯波動。

為探究橋梁位移變化對輪軌力功率譜平穩特性的影響，降低橋梁抗彎彈性模量，增大橋梁橫、豎向變形。撓跨比取L/1500 與L/900，基于不同的撓跨比求動車首輪對輪軌力功率譜。由圖6～圖8 對比可知：橫向力與豎向力功率譜幅值幾乎無變化；增大橋梁變形對影響輪軌力功率譜值影響很小，橋梁位移變化對輪軌力功率譜平穩特性影響不顯著，輪軌力功率譜平穩特性主要受軌道不平順控制。

圖6 動車首輪對輪軌力功率譜Fig. 6 Power spectrum of wheel rail force on first wheel of motor car

圖8 橋梁撓跨比L/900 時的輪軌力功率譜Fig. 8 Wheel-rail force power spectrum of bridge with deflection span ratio of L/900

3.2 脫軌系數極值概率統計

車-橋系統動力響應安全檢驗通常采用3 倍標準差方法，此方法的前提條件是響應分布服從正態分布。脫軌系數服從柯西分布，因此采用3 倍標準差法則顯然不合適。通過前節分析可知，脫軌系數也服從極值分布。本文通過輪軌力功率譜得到大量輪軌力時域樣本，再結合輪軌力均值求得每次脫軌系數時域樣本最大值，進而求得脫軌系數的極值分布。

采用諧波疊加法對輪軌力功率譜進行轉換，得到10 000 組輪軌力時域樣本。求得10 000 個脫軌系數最大值并進行極值分布檢驗。由圖9 可以看出，動車、拖車首輪對脫軌系數最大值可以很好的服從極值I 型分布。置信度99.73%時，動車首輪對脫軌系數的極值為0.281；拖車首輪對脫軌系數的極值為0.385。

圖9 首輪對脫軌系數極值分布Fig. 9 Extreme value distribution of derailment coefficient of first wheel pair

為了進一步驗證極值I 型分布能很好的適用于脫軌系數最大值概率統計，抽取10 000 個樣本中的一部分做極值I 型分布檢驗，從圖10、圖11 可以看出，極值I 型分布能夠很好的適用與脫軌系數最大值概率統計。

圖10 動車首輪對脫軌系數極值分布Fig. 10 Extreme value distribution of derailment coefficient of first wheel pair of motor car

圖11 拖車首輪對脫軌系數極值分布Fig. 11 Extreme value distribution of derailment coefficient of first wheel pair of trail car

采用Monte Carlo 法驗證，其計算樣本數為1000。取每次計算結果的最大值，與極值I 型分布函數的99.73%置信度所對應極值對比。由圖12可知，動車首輪對與拖車首輪對1000 次時程樣本脫軌系數最大值各有2 次超過相應99.73%置信度對應極值。因此運用極值I 型分布求極大值是合理有效的。

圖12 Monte Carlo 法與極值法所得脫軌系數最大值Fig. 12 Maximum value of derailment coefficient obtained by Monte Carlo method and extreme value method

3.3 基于不同車速的脫軌系數極值分布檢驗

列車運行速度對其脫軌系數影響顯著。當波長范圍一定時，速度越大，計算頻率越高，橫、豎向輪軌力越大。因此研究不同速度對脫軌系數的影響以及對極值I 型分布準確性的影響尤為重要。為研究車速變化對脫軌系數的影響，結合實際情況，選取270 km/h、300 km/h、330 km/h，圖13 與圖14 分別給出了極值I 型分布擬合與置信度99.73%的極值與100 次脫軌系數時程樣本最大值關系圖。

對不同車速下脫軌系數分布檢驗可以看出，列車速度越大，置信度99.73%的脫軌系數值越大，且對應的極值分布也越加發散，從圖13(c)可以看出，當速度為330 km/h 時，置信度99.73%的脫軌系數極值為0.438。雖然發散程度存在差異性，但是脫軌系數分布均能很好的擬合于極值I 型分布。

圖13 脫軌系數極值分布Fig. 13 Extreme value distribution of derailment coefficient

由圖14 可看出，對于不同速度，100 次脫軌系數時程樣本最大值均小于置信度99.73%的極值。因此極值I 型分布置信度99.73%所對應的極值可以較好的控制住不同速度下的脫軌系數最大值。

圖14 Monte Carlo 法與極值法所得脫軌系數最大值Fig. 14 Maximum value of derailment coefficient obtained by Monte Carlo method and extreme value method

4 結論

本文采用虛擬激勵法求得輪軌力非平穩功率譜，結合調制函數與諧波疊加法得到輪軌力時域樣本，通過極值理論求得脫軌系數最大值的概率分布，并與Monte Carlo 法對比驗證。算例結算表明本文方法的正確性，并得到以下結論；

(1)本文計算方法能合理地對不同速度下列車脫軌系數進行統計分析，采用極值I 型分布能夠很好的適用于脫軌系數最大值的概率統計。

(2)橋梁變形對輪軌力功率譜平穩特性的影響很小，因此，對于高速鐵路簡支橋梁而言，輪軌力功率譜非平穩特性不顯著。