車輛-軌道非線性耦合動力學(xué)的精細(xì)積分法及其應(yīng)用

張 健，譚述君，2，吳昌華

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116024;2.大連理工大學(xué) 航空航天學(xué)院，大連 116024)

移動車輛與軌道相互作用的數(shù)值模擬可高效準(zhǔn)確地獲得二者的動態(tài)響應(yīng)，為車輛和軌道設(shè)計提供理論依據(jù)。車輛與軌道相互作用時，車輛系統(tǒng)、軌道系統(tǒng)及輪軌關(guān)系存在非線性特性，這時需要在時域下通過逐步積分求解車輛-軌道耦合系統(tǒng)的運(yùn)動方程組。

目前計算車輛-軌道耦合系統(tǒng)的運(yùn)動方程的逐步積分方法主要有:Newmark 法［1-2］，Runge-Kutta 法［3］，Admas法［4］，翟婉明的顯式積分法［5-6］和精細(xì)積分法［7］。文獻(xiàn)［1-6］采用逐步積分方法求解，移動車輛與軌道相互作用時，在每個積分步內(nèi)載荷位置與大小是固定不變的，經(jīng)歷一個積分步長后，突然移到下一個積分步的作用位置。這種計算方法未考慮移動載荷在結(jié)構(gòu)上的連續(xù)移動，載荷中的高頻分量很難準(zhǔn)確獲得［8］。然而利用鐘萬勰提出的精細(xì)積分方法可求解移動載荷在每個積分步內(nèi)連續(xù)變化情況［8-12］。但文獻(xiàn)［8-9，11-12］忽略質(zhì)量慣性作用對結(jié)構(gòu)的影響。全玉云［7］采用精細(xì)積分方法分析車輛-軌道耦合系統(tǒng)動態(tài)響應(yīng);將車輛和軌道的運(yùn)動方程耦合到整體方程中，這導(dǎo)致車輛在軌道運(yùn)行時系統(tǒng)矩陣有時變性。當(dāng)系統(tǒng)矩陣發(fā)生改變時，相應(yīng)的指數(shù)矩陣也要重新計算。張志超［10］提出利用精細(xì)積分方法計算移動質(zhì)量與橋梁相互作用的迭代算法。而其假設(shè)移動質(zhì)量與橋梁始終保持接觸而不發(fā)生分開，無法模擬質(zhì)量與橋梁瞬間脫開現(xiàn)象。文獻(xiàn)［7-12］在求非齊次項(xiàng)的Duhamel積分時需要求系統(tǒng)矩陣的逆陣。譚述君提出非齊次項(xiàng)的Duhamel積分的精細(xì)積分算法［13］，避免計算系統(tǒng)矩陣求逆。

本文吸收他人的研究工作優(yōu)點(diǎn)，提出利用精細(xì)積分求解移動質(zhì)量與結(jié)構(gòu)相互作用的新的方法。采用多項(xiàng)式與指數(shù)矩陣乘法的積分精細(xì)積分算法計算非齊次項(xiàng)的Duhamel積分，避免系統(tǒng)矩陣求逆;利用Lagrange插值多項(xiàng)式獲得載荷項(xiàng)外插公式，避免迭代求解移動質(zhì)量與結(jié)構(gòu)線性/非線性相互作用;分別建立車輛和軌道各個部件的運(yùn)動方程，通過部件之間相互作用關(guān)系把各個部件連接起來，避免當(dāng)車輛在軌道移動時重復(fù)計算各個部件的指數(shù)矩陣。通過有限元形函數(shù)分解移動載荷項(xiàng)，獲得移動載荷在積分步內(nèi)連續(xù)變化的情況。本文建立移動車輛-軌道耦合垂向動力學(xué)模型。模型中考慮車輛與軌道非線性輪軌關(guān)系和軌枕懸空后與道床再次接觸情況。分析鋼軌存在焊接接頭不平順時，軌枕正常支撐狀態(tài)和軌枕懸空對輪軌力和軌道部件之間相互作用力的影響。

1 精細(xì)積分

系統(tǒng)非線性振動方程為:

把式(1)轉(zhuǎn)換到狀態(tài)空間下描述:

其中:

I為單位陣。

利用 Duhamel積分求解式(2)得［16］:

其中T=eHη，η為積分步長。

根據(jù)鐘萬勰［17］提出指數(shù)矩陣的求解方法可以得到T的高精度數(shù)值解。

1.1 非移動載荷項(xiàng)處理方法

由于式(1)中的外力項(xiàng)與系統(tǒng)位移有關(guān)，為避免精細(xì)積分方法迭代求解式(1)，利用Lagrange插值多項(xiàng)式近似獲得外力項(xiàng)。

利用的Lagrange插值多項(xiàng)式獲得在第k時刻到k+1時刻之間的外力項(xiàng)外插公式:

其中:

fk，fk-1和fk-2是第k，k-1和k-2時刻的外力項(xiàng)值。

把式(4)代入式(3)可得:

其中:

利用四階Runge-Kutta法解決精細(xì)積分法起步問題。

1.2 移動載荷項(xiàng)處理方法

當(dāng)載荷的作用位置隨時間變化時，利用形函數(shù)獲得移動載荷在每個積分步內(nèi)連續(xù)變化情況［8］。

圖1中，單元長度為L。在tk時刻，載荷作用位置在A處，經(jīng)過積分步長η后，作用點(diǎn)移到B處，A和B在同一單元內(nèi)。當(dāng)載荷從A運(yùn)行到C處經(jīng)歷τ時間，τ∈［0，η］。A和B距節(jié)點(diǎn)I的距離分別為l1和l2。

圖1 移動載荷作用位置Fig.1 Location of moving load

移動載荷P通過形函數(shù)分解獲得節(jié)點(diǎn)載荷向量:

其中:E為節(jié)點(diǎn)載荷指示矩陣，E隨載荷作用位置而改變。

n為形函數(shù)組成的向量:

把式(7)代入式(9)得:

其中:bij(i=1～4，j=0～3)的具體表達(dá)式見附錄A。

將式(4)和式(10)代入式(8)可得:

其中:ei=E{d1i，d2i，d3i，d4i}T(i=0 ～ 5)

dij(i=1～4，j=0～5)的具體表達(dá)式見附錄B。

把式(11)代入式(3)可得:

其中:T0，T1和T2具體表達(dá)式見式(6)。

如果共有N個移動載荷同時作用在結(jié)構(gòu)上，把每個移動載荷項(xiàng)按上述方法計算非齊次項(xiàng)的Duhamel積分，最后進(jìn)行累加。

當(dāng)積分步長跨越一個單元或多個單元時，計算方法與積分步長不跨越一個單元時類似。假設(shè)積分步起始點(diǎn)和終點(diǎn)構(gòu)成虛擬單元，利用形函數(shù)把移動載荷分配到虛擬單元的節(jié)點(diǎn)上，然后再次利用形函數(shù)把虛擬節(jié)點(diǎn)上的載荷分配到實(shí)際節(jié)點(diǎn)上。具體推導(dǎo)公式見文獻(xiàn)［12］。

1.3 非齊次項(xiàng)的Duhamel積分的算法

鐘萬勰［17］給出式(6)的表達(dá)式:

但是當(dāng)剛度矩陣為奇異矩陣時，無法求解系統(tǒng)矩陣H逆陣。然而根據(jù)譚述君［13］提出的多項(xiàng)式與指數(shù)矩陣相乘的積分的精細(xì)積分算法可以精確獲得T0、T1和T2的數(shù)值解，這樣就避免系統(tǒng)矩陣H的求逆。

構(gòu)造非齊次項(xiàng)的Duhamel積分的加法定理:

利用二項(xiàng)式展開定理式(16)變?yōu)?

將積分步長分為2N份:

這時H(ξ-τ)是非常小的量，可以利用Taylor展開求 eH(ξ-τ)，并取展開項(xiàng)的前 5 項(xiàng):

通過執(zhí)行N次加法定理可以得到多項(xiàng)式與指數(shù)矩陣乘積的積分的精確數(shù)值解。

2 車輛－軌道非線性耦合動力學(xué)模型

車輛-軌道非線性耦合動力學(xué)模型主要包括車輛模型、軌道模型和輪軌關(guān)系。僅考慮車輛和軌道的垂向振動，假設(shè)車輛和軌道結(jié)構(gòu)關(guān)于軌道縱向中心線對稱，因此只建立單側(cè)垂向半節(jié)車輛-軌道耦合動力學(xué)模型，如圖2所示。

圖2 垂向車輛-軌道耦合動力學(xué)模型Fig.2 Vertical vehicle-track nonlinear coupling dynamic model

鋼軌模型采用離散支撐的兩端簡支Euler梁模擬，其余部件簡化成質(zhì)量塊。車輛模型共包括5個自由度:車體沉浮、轉(zhuǎn)向架沉浮與點(diǎn)頭和輪對沉浮。鋼軌考慮橫向和彎曲運(yùn)動。軌枕和道床只考慮垂向運(yùn)動。利用并聯(lián)的彈簧阻尼單元模擬一系和二系懸掛裝置的剛度和阻尼特性、鋼軌與軌枕、軌枕與道床塊、相鄰道床塊之間和道床塊與地基的相互作用。忽略地基的質(zhì)量。

由于沒有考慮鋼軌的阻尼，因此從方程(1)變成鋼軌的運(yùn)動方程為:

其中:Mr和Kr分別為鋼軌質(zhì)量陣和剛度陣;xr和分別為鋼軌的位移和加速度向量;fwr為移動輪軌載荷向量;fr為軌枕對鋼軌作用力向量。

除了鋼軌采用彈性體模擬外，車輛-軌道耦合系統(tǒng)中其他部件采用剛體模擬。

不考慮剛體的阻尼和剛度對剛體運(yùn)動影響，從方程(1)轉(zhuǎn)變成單自由度剛體運(yùn)動方程:

例如轉(zhuǎn)向架的運(yùn)動方程為:

其中:mb為轉(zhuǎn)向架質(zhì)量，為轉(zhuǎn)向架加速度，F(xiàn)c為車體對轉(zhuǎn)向架的作用力，F(xiàn)w為輪對對轉(zhuǎn)向架的作用力。

輪軌法向力由赫茲非線性彈性接觸理論確定［14］:

其中:Fwr為輪軌力，G為輪軌接觸常數(shù)，xw，xr和xirr分別為車輪位移、鋼軌位移和軌道不平順。當(dāng)輪軌法向力為零時，表示車輪與鋼軌脫開。

懸空軌枕與道床采用非線性彈簧單元模擬［15］，具體表達(dá)式為:

其中:Fsb為軌枕與道床相互作用力，kb為道床剛度，xs、xb和dgap分別為軌枕位移、道床位移和軌枕與道床之間間隙。當(dāng)軌枕與道床之間有相互作用力，這時反映軌枕與道床再次接觸情況。

采用余弦函數(shù)模擬焊接接頭不平順［14］:

其中:V為車輛運(yùn)行速度，a為不平順波深，L為不平順波長。

鋼軌運(yùn)動方程的載荷向量包括移動載荷(輪軌力)和固定載荷(軌枕對鋼軌的作用載荷)。移動載荷和固定載荷都是由兩個物體位移決定。例如輪軌力是由車輪與鋼軌位移決定的。然而車輪和鋼軌的位移是由輪軌力決定的。所以在每個時間積分步內(nèi)，輪軌力和車輪與鋼軌的位移是相互耦合的。利用1.1和1.2節(jié)的方法分別可以解除固定載荷和移動載荷與位移的耦合關(guān)系。因?yàn)槔们皫讉€積分步的載荷信息預(yù)測當(dāng)前步的載荷信息，再利用當(dāng)前步的載荷求解車輪和鋼軌的位移。解除載荷與位移之間耦合關(guān)系的優(yōu)點(diǎn)是不需迭代求解鋼軌的運(yùn)動方程。其它部件的運(yùn)動方程處理方法與鋼軌類似。

把方程(22)轉(zhuǎn)換成狀態(tài)空間下的方程(2)形式，其中H矩陣為:

這時無法求解H-1，即無法利用傳統(tǒng)精細(xì)積分求解式(14)中的T0、T1和T2的值。這時需要采用1.3節(jié)提供的方法獲得T0、T1和T2的值。

因此采用顯示精細(xì)積分方法可以無需迭代和求解H-1而計算運(yùn)動方程(21)和(22)。

3 焊接接頭和軌枕懸空對車輛－軌道耦合系統(tǒng)的影響

車輛和軌道具體參數(shù)見文獻(xiàn)［5］。軌道模型長度24 m。焊接接頭屬于上凸型，位于軌道中央的軌枕正上方，波長=0.1 m，波深=0.000 3 m。設(shè)焊接接頭下方軌枕發(fā)生懸空，軌枕與道床之間的間隙為0.001 m。車輛運(yùn)行速度為160 km/h。

圖3為軌枕正常支撐和無軌道不平順時采用精細(xì)積分與Runge-Kutta法的輪軌力計算結(jié)果。兩種結(jié)果非常吻合，這表明本文精細(xì)積分方法的正確性。

圖4是軌枕正常支撐和軌枕發(fā)生懸空時，圖2中右側(cè)輪軌力隨時間變化情況。在車輪未進(jìn)入焊接接頭不平順和軌枕懸空的影響區(qū)內(nèi)時，軌枕正常支撐和軌枕懸空時輪軌力相同(圖4未具體給出)。在車輪未進(jìn)入焊接接頭不平順影響區(qū)，但進(jìn)入軌枕懸空的影響區(qū)時，兩種情況的輪軌力略有差別，主要是軌枕懸空后，軌道的支撐剛度變低。由于只有一根軌枕發(fā)生懸空，它對軌道支撐剛度影響范圍很小，所以這時兩種輪軌力差別不大。當(dāng)車輪進(jìn)入焊接接頭不平順的影響區(qū)內(nèi)時，軌枕懸空時輪軌力的第一個波峰略大于軌枕正常支撐時的輪軌力。這是因?yàn)榈谝粋€輪軌力波峰是主要由焊接接頭不平順造成的，懸空軌枕對輪軌力影響較小，所以兩種情況的輪軌力相差不大。輪軌力的最大值發(fā)生在焊接接頭不平順深度逐漸增加時不平順的最大斜率附近，這時不平順深度隨時間變化非常快。經(jīng)過斜率最大值后，不平順深度增加速度漸漸降低，輪軌力也隨著降低。當(dāng)不平順深度達(dá)到最高值時，斜率為零，而輪軌力基本恢復(fù)到無焊接接頭不平順時情況。不平順深度漸漸降低和較大輪軌力，造成車輪與軌道發(fā)生脫離現(xiàn)象。兩種情況下車輪和鋼軌脫離時間幾乎相同，這是因?yàn)閮煞N情況的第一個輪軌力波峰相差不大。由于車輛重力影響，車輪再次與鋼軌發(fā)生接觸，輪軌力產(chǎn)生第二個波峰。這個輪軌力波峰值小于第一個輪軌力波峰值，軌枕懸空時的輪軌力明顯大于軌枕正常支撐時的輪軌力。因?yàn)檫@時車輪已經(jīng)過焊接接頭不平順影響區(qū)域，但是還在軌枕懸空區(qū)影響區(qū)域內(nèi)，輪軌力差異由軌道的支持剛度影響。由于輪軌力第二個波峰持續(xù)時間較長，即沖擊力從車輪與鋼軌接觸界面向車輛和軌道傳播時間長，會對遠(yuǎn)離車輪與鋼軌接觸界面的車輛和軌道的構(gòu)件造成一定損傷。車輪經(jīng)過軌枕懸空影響區(qū)后，由于車輛和軌道阻尼的作用，輪軌力漸恢復(fù)到未經(jīng)歷軌道不平順情況。圖2中左側(cè)輪軌力與右側(cè)的相似，這里就不再介紹。

圖3 精細(xì)積分與Runge-Kutta法的輪軌力計算結(jié)果對比Fig.3 Comparsion of wheel/rail force calculated by precise integration method and Runge-Kutta

圖4 輪軌力時間歷程Fig.4 Wheel/rail force time history

圖5 輪軌力、鋼軌與軌枕相互作用力和軌枕與道床相互作用力的最大值Fig.5 The maximum wheel/rail force，interaction force of the rail and the sleeper and interaction force of the sleeper and the ballast

在右側(cè)車輪作用下，輪軌力、鋼軌與軌枕相互作用力和軌枕與道床相互作用力的最大值發(fā)生的運(yùn)行時間，如圖5所示。在軌枕正常支撐時，當(dāng)輪軌力從車輪與鋼軌接觸界面向下傳時，由于阻尼吸收能量，它們數(shù)值依次減小，而且有一定時間差。當(dāng)軌枕與道床發(fā)生再次接觸，它們之間的相互作用力比較小。這是因?yàn)閼铱哲壵淼膬蓚?cè)軌枕承擔(dān)大部分來自輪軌的作用力。

圖6是右側(cè)車輪通過焊接接頭不平順時，焊接接頭下方軌枕及其左右各一根軌枕與道床相互作用力比例因子圖。比例因子是軌枕正常支撐時或軌枕懸空時的軌枕與道床相互作用力與軌道無任何缺陷時的比值。軌道存在接頭不平順時，焊接接頭正下方的軌枕與道床相互作用力為軌道正常時的1.6倍。這就容易造成道床下沉，最終導(dǎo)致道床失去對軌枕的支撐能力，出現(xiàn)軌枕懸空情況。當(dāng)軌枕懸空出現(xiàn)后，鄰近懸空區(qū)的軌枕與道床作用力明顯增大，這就容易造成軌枕懸空區(qū)進(jìn)一步擴(kuò)大。

圖6 軌枕與道床相互作用力Fig.6 Interaction force of the sleeper and the ballast

4 結(jié)論

利用本文方法求解移動質(zhì)量與結(jié)構(gòu)相互作用的特點(diǎn)是不需要求解系統(tǒng)矩陣逆陣，不需要迭代求解考慮移動質(zhì)量和結(jié)構(gòu)相互作用的非線性和結(jié)構(gòu)本身的非線性的運(yùn)動方程，系統(tǒng)矩陣具有時不變性，而且可以考慮移動質(zhì)量在每個積分步內(nèi)連續(xù)變化情況。

分析車輛在軌道運(yùn)行時，輪軌力和軌道部件之間的相互作用力的時間歷程。計算結(jié)果表明，在焊接接頭不平順的影響區(qū)內(nèi)，一根懸空軌枕對輪軌力影響不大。軌枕懸空主要影響第二個輪軌力峰值。焊接接頭不平順引起較大的輪軌力易造成道床下沉，出現(xiàn)軌枕懸空現(xiàn)象。出現(xiàn)一根軌枕懸空后，易造成鄰近軌枕也發(fā)生懸空現(xiàn)象。

本文方法可以進(jìn)一步應(yīng)用到車輛與無碴軌道相互作用、車輛與橋梁相互作用和汽車與路面相互作用等移動質(zhì)量與結(jié)構(gòu)相互作用的領(lǐng)域。

［1］Lei X，Noda N A.Analyses of dynamic response of vehicle and track coupling system with random irregularity of track vertical profile［J］.Journal of Sound and Vibration，2002，258(1):147-165.

［2］Zhang Y W，Lin J H，Zhao Y，et al.Symplectic random vibration analysis of a vehicle moving on an infinitely long periodic track［J］.Journal of Sound and Vibration，2010，329(21):4440-4454.

［3］ Li Z G，Wu T X.On vehicle/track impact at connection between a floating slab and ballasted track and floating slab track's effectiveness of force isolation［J］.Vehicle System Dynamics，2009，47(5):513-531.

［4］Nielsen J C O，Igeland A.Vertical dynamic interaction between train and track-influence ofwheeland track imperfections［J］.Journal of Sound and Vibration，1995，187(5):825-839.

［5］ Jin X S，Wen Z F，Wang K Y，et al.Three-dimensional train-track model for study of rail corrugation［J］.Journal of Sound and Vibration，2006，293(3-5):830-855.

［6］Zhai W M，Wang K Y，Cai C B.Fundamentals of vehicletrack coupled dynamics［J］.Vehicle System Dynamics，2009，47(11):1349-1376.

［7］全玉云.機(jī)車車輛/軌道系統(tǒng)垂向耦合動力學(xué)有限元分析的研究［D］.北京:鐵道部科學(xué)研究院，2000.

［8］林家浩，張守云，呂 峰，等.移動簡諧載荷作用下橋梁響應(yīng)的高效計算［J］.計算力學(xué)學(xué)報，2006，23(4):385-390.

［9］余 華，吳定俊，項(xiàng)海帆.移動荷載過橋的精細(xì)計算［J］.振動與沖擊，2009，28(5):17-21.

［10］張志超，林家浩.移動質(zhì)量作用下橋梁響應(yīng)的精細(xì)積分［J］.大連理工大學(xué)學(xué)報，2006，46(S1):26-32.

［11］ Lu F，Kennedy D，Williams F W，et al.Non-stationary random vibration of fe structures subjected to moving loads［J］.Shock and Vibration，2009，16(3):291-305.

［12］呂 峰，林家浩，張亞輝.移動隨機(jī)荷載作用下橋梁振動分析［J］.振動與沖擊，2008，27(12):73-78.

［13］譚述君，鐘萬勰.非齊次動力方程duhamel項(xiàng)的精細(xì)積分［J］.力學(xué)學(xué)報，2007，39(3):374-381.

［14］翟婉明.車輛-軌道耦合動力學(xué)(第三版)［M］.北京:科學(xué)出版社，2007.

［15］Bezin Y，Iwnicki S D，Cavalletti M，et al.An investigation of sleeper voids using a flexible track model integrated with railway multi-body dynamics［J］. Proceedings ofthe Institution of Mechanical Engineers Part F-Journal of Rail and Rapid Transit，2009，223(6):597-607.

［16］鐘萬勰.應(yīng)用力學(xué)的辛數(shù)學(xué)方法［M］.北京:高等教育出版社，2006.

［17］Zhong W X，Williams F W.A precise time-step integration method［J］.Proceedings of the Institution of Mechanical Engineers PartC-JournalofMechanicalEngineering Science，1994，208(6):427-430.

附錄A

bij系數(shù)表達(dá)式:

附錄B

dij系數(shù)表達(dá)式

同理可得其它值。