劉鵬飛 *, 楊紹普 *, 劉永強 ** 顧曉輝 *, 劉澤潮 *,

* (石家莊鐵道大學省部共建交通工程結構力學行為與系統(tǒng)安全國家重點實驗室,石家莊 050043)

? (石家莊鐵道大學河北省交通工程結構力學行為演變與控制重點實驗室,石家莊 050043)

** (石家莊鐵道大學機械工程學院,石家莊 050043)

引言

輪對是軌道交通車輛最為重要的走行部件,發(fā)揮著承載、牽引、制動及導向等重要功能,但同時屬于簧下質量,在長期使用過程中,車輪踏面及軌面不可避免會出現剝離掉塊、擦傷、波磨及車輪多邊形等典型的短波不平順激擾,容易激發(fā)輪對的中高頻結構振動,甚至輪對結構與軌道結構的自振頻率相吻合而演化為系統(tǒng)的耦合共振,由此惡化了列車的運行環(huán)境并加劇輪軌缺陷的進一步發(fā)展.因此在車輛系統(tǒng)中,中高頻域動力學研究中考慮輪對結構振動十分必要,對于延緩和消除相關問題的產生具有重要意義.

近年來,從輪對柔性建模及動態(tài)仿真問題,國內外相關學者開展了大量研究.主要的研究方法是,基于有限元軟件進行輪對彈性建模,提取振動模態(tài),之后與動力學軟件進行聯合仿真,進而實現非線性振動的求解.例如,文獻[1-12]主要建立了輪對的有限元模型,進一步與動力學仿真軟件結合,從而開展中高頻輪對合理建模問題[1]、輪軌磨耗問題[2]、多邊形激擾問題[3-4]、車輪扁疤沖擊振動問題[5-6]、車輛運動穩(wěn)定性問題[7]、鋼軌波磨激擾問題[8-9]、軸承壽命問題[10]、滾動噪聲問題[11]、輪軌摩擦自激振動與輪對彎曲振動耦合問題[12]等方面的研究.可以看出,相關的建模方法和仿真技術也是當前剛柔耦合動力學的主要研究手段.此外,國內外學者也在積極探索輪對彈性振動建模的其他方法.如Guiral 等[13]在非慣性車輛運動參照系的框架內推導出了一種柔性輪對的計算公式.Baeza 等[14]給出了三種基于歐拉方法的柔性轉動固體動力響應計算公式,適用于固體與非旋轉結構相互作用的研究.徐寧等[15]利用假設模態(tài)法,得到帶有集中質量及轉動慣量的彈性車軸振型函數,從而建立考慮輪對彈性振型的車軸模型.文獻[16]利用歐拉坐標系中的拉格朗日方程和模態(tài)疊加法,建立和求解輪對彈性振動的運動方程.崔瀟等[17]采用歐拉坐標系下旋轉效應柔性輪對模型進行剛柔耦合車輛-軌道系統(tǒng)多體動力學建模,楊云帆等[18]采用歐拉梁模擬輪軸彈性振動,建立了考慮輪對柔性的直線電機地鐵車輛動力學模型.

總體而言,輪對彈性振動的建模方法較為成熟,有限元方法最為常用,通過模態(tài)疊加原理實現隨機激擾下輪對結構振動與車輛振動的求解,結合有限元和動力學軟件是一種有效的途徑.但無論是基于有限元法還是直接推導輪對的彈性振動方程,結構改變后往往需重新建模和推導,較難實現模塊化,計算機編程也略為復雜.近年來,傳遞矩陣法已突破線性振動的限制,通過時間離散實現了轉子系統(tǒng)、武器系統(tǒng)等非線性動力學及彈性振動問題的求解[19-22],該方法通過可將振動系統(tǒng)分解為一定數量的離散單元,各單元均可獨立編制傳遞矩陣,最后再通過單元邊界處的力和運動狀態(tài)向量集成為整體振動系統(tǒng),無需推導系統(tǒng)總體動力學方程,因而具有求解規(guī)模少、模塊化組合、易于編程等優(yōu)點,但相關方法在軌道交通車輛動力學研究中還鮮有報道.

因此,本文以機車車輛單軸滾振試驗臺為研究對象,引入分布質量彈性輪軸的離散時間傳遞矩陣、車輪傳遞矩陣等,進一步建立考慮輪對彈性振動的試驗臺剛柔耦合動力學分析模型,綜合新型顯式積分法、Newmark 隱式積分法和Riccati 法對系統(tǒng)振動進行數值求解,最后通過單軸滾振試驗臺的測試結果對仿真模型進行驗證,以期為考慮結構部件彈性振動的車輛-軌道耦合系統(tǒng)動力學研究提供新的建模方法和研究思路.

1 分布質量彈性輪軸離散時間傳遞矩陣

用傳遞矩陣法對輪對建模時,可將輪對分為分布質量彈性輪軸和集總質量的車輪等剛體部件,通過狀態(tài)向量與傳遞矩陣的關系實現截面力和運動量的傳遞計算.狀態(tài)向量用于描述部件端部的剪力、彎矩、位移及加速度等力和運動量.文獻[19]中系統(tǒng)建立了以力、位移為狀態(tài)變量的傳遞矩陣法.文獻[21]針對無質量彈性梁-滑動軸承-圓盤系統(tǒng)進行了研究,同時指出,研究軸系結構瞬態(tài)響應時,由于彈性模量大、時間步長小,可能出現病態(tài)矩陣,采用力、加速度的狀態(tài)向量可盡力避免此問題.鑒于此,首先設定輪軸狀態(tài)變量W為

式中,M為截面力矩,Q為剪力,為截面轉角加速度,為截面的垂向振動加速度.

建立分布質量彈性輪軸離散時間傳遞矩陣的本質在于確定彈性軸輸入端和輸出端狀態(tài)向量的力學傳遞關系.對于動力輪對結構(圖1),從參振的角度可大體分為軸箱、車輪、齒輪箱與輪軸四個部分,其中輪軸作為相對的細長結構,貢獻了主要的彈性振動能量.彈性輪軸可視為彈性梁結構,假設其第i段長度為li,梁的線密度為ρ,在車輛坐標系中定義左側為輸入端,剪力QL,彎矩ML,定義右側為輸出端剪力QR、彎矩MR,則其受力狀態(tài)如圖2 所示,對應的左右側位移為ZL和ZR,g為重力加速度,ρg為梁上分布質量的重力分布.

圖1 輪對結構Fig.1 Wheelset structure

圖2 分布質量彈性軸的受力與變形Fig.2 Force and deformation of elastic axle with distributing mass

設梁的彈性振動位移用模態(tài)疊加法求解,引入正則振型函數φ(x),正則振型坐標q(t),則梁在距其左端面x處、t時刻的彈性位移可表示為

根據其剪力和彎矩平衡關系,并考慮慣性力,可建立梁單元的受力平衡方程式(3)和式(4),對于輸入端到輸出端的垂向位移及轉角變形,同樣根據剪力、彎矩及分布載荷條件下的梁彈性變形進行疊加可得到變形方程式(5)和式(6)[22]

式中,EI為梁的抗彎剛度.式(3)～ 式(6)即為輸入端到輸出端的分布質量彈性梁的力和位移傳遞關系.為了進行數值積分求解同時保證彈性梁振動求解的收斂性和穩(wěn)定性,引入Newmark-β隱式積分法進行上述四式的改造[19].Newmark-β法積分格式如下

基于式(7)～ 式(10),將式(3)～ 式(6)中的位移項表示為當前時刻加速度及上一時刻加速度、位移及振動速度的顯式表達,則改造后的位移公式可進一步表達為式(11)和式(12)

上式中還存在積分項,進行傳遞矩陣推導的關鍵在于消除模態(tài)函數積分項中的正則加速度坐標.以式(3)中積分項為例,若取前4 階模態(tài),結合式(2),可知

根據振動力學,梁端位移和轉角的加速度滿足如下關系

于是可將轉角、位移的加速度統(tǒng)一用矩陣表示為式(15),式中U為可用模態(tài)函數求解出的系數矩陣

令

則關于振型函數的定積分項可表示為

對比式(15)和式(16),加速度狀態(tài)向量與積分項的矩陣表達都是關于正則坐標加速度的.那么,如果通過加速度狀態(tài)向量的線性疊加組合來表示該積分項,就可消除原方程積分項中關于時間的正則振型坐標加速度.為此,引入待定系數a1～a4,則積分項可進一步表達為

式中,左端和右端項可約掉正則振型坐標加速度項,左右側再經過轉置,可得

最終通過矩陣運算,可得待定系數a1～a4的具體數值

式(18)中,方程的左側是正則振型函數關于軸長度的積分,在已知振型函數后,可直接求得結果,同樣轉置矩陣UT也可通過振型函數獲得,于是關于待定系數a1～a4的求解就轉化成了一階常系數線性方程組的求解問題.只要振型假設合理,矩陣UT可逆,便可求出待定系數的具體量值,也即表明通過軸左、右端的狀態(tài)向量間接表達了積分項.對于式(18)也可通過高斯消元法快速求解.

同理,對于式(4)、式(11)和式(12)中的包含振型函數的積分項,采用類似方法也可得到對應的組合待定系數bi,ci和di(i=1～ 4).最終,將方程式(3)～式(4)、式(11)～式(12),統(tǒng)一表達為如下矩陣形式

則,左右端的傳遞關系可轉變?yōu)?/p>

式中,H和D′矩陣為軸段的場傳遞矩陣,A1和A2及D為推導出的中間矩陣,且滿足下式

2 車輪、齒輪箱及彈簧阻尼系統(tǒng)傳遞矩陣

對于車輪,簡化為具有垂向、側滾位移的剛體處理,其左、右側受力狀態(tài)及狀態(tài)向量如圖3 所示.剛體質心兩側的轉角、位移均相同,基于達朗貝爾原理,根據垂向力平衡和轉矩平衡關系,不難推導出其振動微分方程,此處不再贅述,直接給出車輪的傳遞矩陣,如式(25)中Hw和Pw.

圖3 車輪受力狀態(tài)Fig.3 Force conditions of wheel

式中,Jw為車輪側滾慣量,mw為車輪質量,Fw為輪軌垂向力,Mw為側滾蠕滑力矩.

對于齒輪箱,將其2/3 質量作為簧下參振質量,與彈性輪軸連接,則式(25)中傳遞矩陣同樣適用于齒輪箱.對于軸箱,僅考慮軸箱的垂向振動,為此式(25)中進一步令Jw和Mw=0 即可.

轉向架系統(tǒng)中,軸承、一系和二系懸掛均可一定程度上簡化為剛度-阻尼力元,本文以軸承為例進行說明.如圖4 所示的彈簧、阻尼系統(tǒng),彈簧剛度為Kb,阻尼系數Cb,左右側的位移不相等,剪力相同、方向相反,若不考慮該系統(tǒng)的彎曲剛度和阻尼,可認為左右側的彎矩相同,則滿足式(26)中的數學關系

圖4 彈簧-阻尼系統(tǒng)受力狀態(tài)Fig.4 Force conditions of spring-damping system

進一步,將式(26)采用矩陣表示,并引入Newmark-β法積分格式對其改造,最終可得到式(27)的表達形式,即為該力元的離散時間傳遞矩陣

3 單軸滾振試驗臺的垂向動力學建模及求解方法

3.1 動力學建模

以單軸滾振試驗臺為例(圖5),建立其動力學模型(圖6).需要指出的是,模型分為兩部分,一為輪對系統(tǒng),基于上文的離散時間傳遞矩陣法進行建模;二為虛擬構架和軌道輪組部分,直接給出動力學方程,并采用新型顯式積分法(翟方法)[23]求解其速度和位移響應.這樣處理的目的是,輪對之外的系統(tǒng)采用顯式積分法求解容易編程實現,且無需聯立求解高階代數方程組,可進一步提高計算效率.建模時,可將輪對彈性體作為子模型,方便后期與車輛-軌道耦合動力學模型、列車-軌道動力學模型[23-26]及其求解方法進行對接,同時也可通過本次分析和試驗,來檢驗復雜剛柔耦合系統(tǒng)中顯、隱式混合積分求解振動的可行性.

圖5 單軸滾振試驗臺系統(tǒng)構成Fig.5 Compositions of single-wheelset rolling and vibration test system

輪對發(fā)揮彈性振動的主要結構是輪軸,可將輪軸采用分布質量彈性梁建模.軸承的支撐剛性和阻尼則簡化為彈簧-阻尼力元.如圖6 所示,本文將輪對-齒輪箱系統(tǒng)從左到右分為12 個部分,3,5,6,8,10 為彈性軸段,采用式(21) 中傳遞矩陣;軸箱體1 和12 采用式(25)傳遞矩陣,但僅考慮垂向振動;軸承力元2 和11 采用式(27)中的傳遞矩陣;車輪4 和9 及齒輪箱7 采用式(25)的傳遞矩陣.

圖6 試驗臺動力學模型Fig.6 Dynamic model of test rig

對于剛性構架和軌道輪,僅考慮其垂向運動,對應的振動微分方程為

由于針對垂向振動,輪軌垂向力簡化為Hertz 非線性接觸彈簧[23],以右輪為例,作用力可表示為

式中,R為車輪名義滾動半徑,Z0(t)為t時刻軌道不平順,ZtR(t)為軌道輪垂向位移.誠然,輪-輪接觸和輪-軌接觸有一定差異,若考慮空間接觸關系,需進一步描述輪-輪接觸斑形狀、蠕滑力的不同,具體可參考文獻[27].

3.2 振動求解方法

下面重點闡述基于Riccati 法來求解輪對系統(tǒng)的振動加速度.基于離散時間傳遞矩陣法建立的各部件傳遞矩陣均具有如式(31) 的統(tǒng)一形式,其中,H11,H12,H21,H22為2 × 2 階方陣,F,E為2 × 1階列陣.將狀態(tài)變量W分為力狀態(tài)變量和加速度狀態(tài)變量兩部分,分別對應f和,如式(32)所示

引入Riccati 變換[21,23],第i截面端部力f和加速度狀態(tài)變量之間滿足如下關系

式中,各單元對應的Pi和Si矩陣可由下式遞推計算

依據系統(tǒng)的邊界條件確定式(33)～ 式(36)的邊界狀態(tài).對于左側軸箱,邊界條件滿足f=0,=0,那么可知P0=0 且S0=0,由式(35)和式(36)依次遞推出各截面的P和S矩陣.進而,由最右端滿足f13=0,=0,則通過下式(37)確定最右端的加速度矢量.最后根據式(33)和式(34)依次獲得各截面的加速度和力狀態(tài)變量

本文模型的求解涉及了3 種重要的計算方法,系統(tǒng)總體求解流程如圖7 所示.具體的求解步驟為:

圖7 動力學模型數值求解流程Fig.7 Numerical solving flow of dynamic model

(1) 獲取t時刻的振動位移、速度,以及輪軌力、軸承力等線性及非線性作用力;

(2)采用新型顯式積分法(翟方法)計算剛性構架、軌道輪t+Δt時刻的位移和速度響應,預估柔性輪對系統(tǒng)振動位移和速度;

(3)基于第(1)步結果和第(2)步預估結果,計算t+Δt時刻一系懸掛力、輪軌力、軸承力;

(4)采用Riccati 方法計算柔性輪對系統(tǒng)振動加速度;

(5) Newmark-β法計算t+Δt時刻柔性輪對系統(tǒng)的振動速度、位移;

(6)檢驗第(2)步預估振動量和Newmark-β法積分求解振動量的差異是否小于允許誤差,若大于設定誤差,進行預估振動量與積分振動量的平均,再返回步驟(3),重新計算懸掛力、輪軌力;重復上述步驟,直到誤差低于設定閾值,即獲得需要的輪軸系統(tǒng)振動速度、位移;

(7)反饋輪對系統(tǒng)振動量及一系懸掛力、輪軌力至構架、軌道輪組成的外部系統(tǒng),通過振動方程直接計算下一時刻構架、軌道輪的加速度;

(8)保存t+Δt時刻系統(tǒng)的位移、速度和加速度,作為下一時刻積分的初始條件.

4 輪對振動特性數值仿真和試驗驗證

針對單軸滾振試驗臺及輪對進行數值仿真和同步試驗.基于假設模態(tài)法,文獻[22]確定振型函數,并代入到前文分布質量彈性軸的傳遞矩陣中.建模涉及的主要參數在表1 中列出,其中軸承的剛度和阻尼可參考文獻[28]的計算結果,軌道輪組下部其支承剛度和阻尼參考文獻[29].至于柔性輪對的自振特性,可直接參照文獻[30]中作者前期給出的計算方法,采用傳遞矩陣亦可方便求解,對于本文所研究的動力輪對,其所考慮的最高階模態(tài)為3 階彎曲模態(tài),對應頻率為560 Hz.

表1 主要建模參數Table 1 Main modeling parameters

4.1 軌道輪切削、打磨與測試

為了檢驗理論模型的正確性,開展了300～400 km/h 速度級輪對激振試驗,軌道輪直徑1.8 m,車輪名義滾動直徑0.915 m.對車輪和軌道輪進行了表面幾何狀態(tài)測試(圖9),并對軌道輪進行了多邊形打磨.測試結果顯示(圖10),輪對有初始的一階不圓順,軌道輪即使通過切削加工也存在較低水平的粗糙度,但幅值不超過0.03 mm.進一步通過局部打磨模擬了軌道輪多邊形和局部凹陷.其中多邊形的波長157 mm,對應于軌道輪36 階多邊形,幅值不超過0.05 mm,局部凹陷波長300 mm、波深0.3 mm,代表單一的諧波不平順,軌道輪切削后初始的軌面粗糙度水平較低.在試驗中,車輪不圓順始終參與系統(tǒng)激振,為此在仿真中將車輪多邊形作周期延拓,并與軌道輪3 種激擾進行疊加,形成復合不平順(圖10(c)),以此作為仿真模型的輪軌接觸界面輸入激勵.

圖9 室內試驗Fig.9 Laboratory test

圖9 室內試驗(續(xù))Fig.9 Laboratory test (continued)

圖10 輪對及軌道輪表面幾何狀態(tài)Fig.10 Geometric states of wheelset and rail roller surfaces

4.2 振動特性的試驗驗證

基于前述離散傳遞矩陣法建立的彈性軸單輪對動力學模型,進行了理論模型的試驗驗證,試驗過程中,測試了不同軌面狀態(tài)下的軸箱振動加速度,信號采集采用壓電式加速度傳感器,采樣頻率不低于5000 Hz,以盡量涵蓋主要的中高頻振動.圖11 給出了初始不平順狀態(tài)下,輪對以400 km/h 速度運行時軸箱加速度測試結果和仿真結果的比較.從頻域上看,仿真結果較好捕捉了軸箱1600 Hz 范圍內的振動主頻峰值,但主要振動能量集中于500 Hz 以內.在時域中,對測試結果進行1600 Hz 和500 Hz 低通濾波,與計算結果進行比較.500 Hz 低通濾波后,測試加速度和仿真加速度幅值總體分別分布在(-3.4～+2.8)g和(-3.6～+3.1)g的范圍內,最大幅值相差0.3g左右,約9%.

圖11 初始粗糙度狀態(tài)下仿真和測試加速度比較Fig.11 Comparisons between simulation and test acceleration under initial roughness

進一步,對多邊形激擾下的軸箱振動響應進行對比驗證,測試和仿真加速度對比如圖12 所示,輪對運行速度400 km/h.同樣,在頻域上看,仿真結果較好反映了實際系統(tǒng)在1510 Hz 范圍內主要振動特征,振動主頻包含了輪對轉頻(39 Hz)、多邊形激擾頻率(707 Hz)及其倍頻等.在時域響應上,對測試結果進行1510 Hz 和1000 Hz 的低通濾波,幅值分布在 ± 8.5g和-6g～ +7.5g的范圍內,仿真加速度幅值總體在-6g～ +7.9g的范圍內,因此在1000 Hz 頻率內,測試和計算結果幅值相差0.4g,約5%.事實上,在500～ 1000 Hz 的頻率范圍內,僅軌道輪多邊形激擾引起的強迫振動較為顯著,其余頻率對應振動能量十分微弱.

圖12 多邊形激擾下仿真和測試加速度比較Fig.12 Comparisons between simulation and test accelerations under polygon excitation

最后,以軌道輪局部凹陷為主要激勵,對仿真模型的計算結果進行了對比驗證,由于局部凹陷的幅值較大,考慮運行安全,在試驗中將輪對走行速度設定在了300 km/h.同時為研究剛、柔建模方法對系統(tǒng)動態(tài)響應帶來的影響,也分析了剛性輪對的振動響應,在時域-頻域對仿真和實測軸箱振動加速度進行了比較,如圖13 所示.對于柔性輪對,仿真結果反映了實際系統(tǒng)在1580 Hz 范圍內主要振動特征,但在500 Hz～ 1580 Hz 的高頻區(qū)段,振動幅值較小,振動能量總體集中于500 Hz 以內.經過1580 Hz 和500 Hz 的低通濾波后,在時域響應上,測試和和仿真結果能清晰識別和捕捉局部凹陷產生的振動激擾成分,加速度幅值范圍均在(-19～+21)g之間,響應規(guī)律較接近,最大幅值相差在1g之內,誤差在5%左右.對于剛性輪對,總體能夠反映250 Hz 以內的振動,但超過該頻率的高頻部分未能充分體現,相較測試結果和柔性輪對計算結果要小,最終在時域響應上,剛性輪對建模的振動加速度幅值分布在 ± 11g范圍內,約為實測結果的60%左右,因此基于離散時間傳遞矩陣法的柔性輪對模型具有更高的求解精度.也進一步表明,采用柔性輪對建模方法能夠更好的描述結構在中高頻域的彈性振動.

圖13 局部凹陷激擾下仿真和測試加速度比較Fig.13 Comparisons between simulation and test accelerations under local-dent excitation

綜上所述,基于離散時間傳遞矩陣法建立的單輪對彈性體模型能夠較好模擬輪軌高頻激振條件下的振動性能,與單軸滾振試驗臺的測試結果對比表明,模型在500 Hz 頻率范圍內的振動能量集中區(qū)總體具有較好的適應性.采用新型顯式積分法、Newmark-β隱式法對系統(tǒng)的混合積分數值求解方案也是可行的,獲得了較好的求解精度.就車輪多邊形磨耗而言,根據以往研究,輪對的彎曲振動是其重要誘因之一[31],所建立的理論模型能夠反映這一要素.但仍需要指出的是,輪對同時存在扭轉、彎曲耦合振動,輪軌法向接觸非線性、蠕滑非線性特性同時會影響輪對的結構振動.誠然,在本文研究的工況下以垂向振動為主,后期在車輛-軌道耦合動力學模型進一步集成該模型時,需要綜合考慮輪對系統(tǒng)的空間彈性振動形態(tài).另外,離散時間傳遞矩陣法可實現軸、盤類等形狀規(guī)則零件的模塊化建模,但如何利用該方法對車輪、齒輪箱等空間復雜曲面及結構等高維問題進行動力學建模仍需深入探索.

5 結論

本文以機車車輛單軸滾振試驗臺為研究對象,推導了彈性輪軸及主要參振部件的離散時間傳遞矩陣,同時考慮上部構架和下部軌道輪組的振動,建立了試驗臺動力學分析模型,綜合新型顯式積分法、Newmark-β隱式積分法及Riccati 方法對系統(tǒng)進行了數值求解.通過打磨軌道輪,開展了高速條件下的輪對運行試驗和仿真測試,主要結論如下.

(1)綜合對比軌道輪微幅不平順、多邊形沖擊和凹陷式諧波激擾下的計算和試驗結果,理論模型總體能夠捕捉及反映輪對結構中高頻的振動特性,特別是在500 Hz 以內的振動能力集中區(qū),時域上的仿真結果響應規(guī)律和測試結果比較接近,三種工況下計算的加速度幅值誤差分別低于9%,5%和5%.

(2)輪對彈性系統(tǒng)采用離散時間傳遞矩陣法建模,采用Riccati 方法和Newmark-β隱式積分法求解,外部構架及軌道輪振動則采用顯式積分求解,上述混合建模和混合積分求解的思路是可行的.

(3)離散時間傳遞矩陣方法建模時無需推導系統(tǒng)總體動力學方程,只要具備了各組成部件的傳遞矩陣,不同模型可調用相應模塊并組裝即可,且涉及的矩陣階次低,可實現快速建模及求解.