周成博,馬道林

(1.西南交通大學(xué)高速鐵路線路工程教育部重點實驗室,成都 610031； 2.西南交通大學(xué)土木工程學(xué)院,成都 610031)

1 研究背景

輪軌黏著規(guī)律是輪軌關(guān)系研究的核心基礎(chǔ)。18世紀(jì)末Carter[1]首先解決了二維彈性體滾動接觸問題。他利用庫倫摩擦定律和Hertz接觸理論成功地計算出接觸斑內(nèi)切向應(yīng)力的大小，劃分出了黏著區(qū)和滑動區(qū)，為輪軌黏著的研究奠定了理論基礎(chǔ)。Vermulen和Johnson[2]則將Carter理論擴展到三維無自旋橢圓接觸。在此基礎(chǔ)上Kalker[3]建立了至今為主應(yīng)用最為廣泛的由線性理論、簡化理論和精確理論組成的蠕滑理論。

以上涉及的都是不考慮輪軌表面粗糙度則在且向上均采用庫倫摩擦模型即摩擦系數(shù)恒定不變的算法。為了更為真實地模擬輪軌接觸，一些學(xué)者考慮使用函數(shù)型摩擦系數(shù)代替庫倫摩擦，如Polach[4]在其模型中引入速度型摩擦函數(shù)以考慮輪軌接觸剪切剛度；張衛(wèi)華等[5]將速度型摩擦系數(shù)引入FASTSIM算法中，常崇義[6]利用ABAQUS軟件建立穩(wěn)態(tài)輪軌滾動接觸有限元模型，其中摩擦系數(shù)隨著輪軌微滑速度變化。

以上方法考慮到了摩擦系數(shù)隨速度的變大，但仍將車輪與鋼軌表面視作光滑表面，導(dǎo)致計算結(jié)果與現(xiàn)場試驗不同。因此一些學(xué)者開始考慮將輪軌表面視作粗糙表面來進(jìn)行接觸計算。孫瓊[7]基于GreenWood的粗糙表面彈性接觸模型，改良后得到粗糙表面的彈-塑性接觸模型，進(jìn)行粗糙輪軌表面的接觸計算；吳兵等[8]則通過數(shù)值模擬的方法計算了水油等第三介質(zhì)存在下的輪軌粗糙表面接觸計算。以上計算仍基于Greenwood[9]對于粗糙度的計算方法，從統(tǒng)計學(xué)的角度出發(fā)，缺乏物理本構(gòu)模型的支持。

本文嘗試從物理學(xué)的角度出發(fā)，旨在研究干燥條件下中尺度的軌道接觸，以較小規(guī)模且隨機的方式理解輪軌摩擦行為。Braun[10]引入中尺度摩擦模型來解釋摩擦的起始，即摩擦從黏著到滑動的過渡。Braun的模型[10]解釋了很多關(guān)于黏滑過渡的實驗結(jié)果，這些結(jié)果是在高于10 000 frames/s的采樣率下用高速攝像機獲得的[11-13]。而輪軌接觸斑內(nèi)滑動區(qū)與黏著區(qū)的存在，則使得與上述論文中類似的方法來模擬輪軌接觸提供了理論基礎(chǔ)，不僅能考慮滑動而且能考慮滾動。Aghababaei[14]對從接觸到兩個粗糙表面之間的粗糙度分離的整個過程的數(shù)值研究與分子動力學(xué)模擬證實了模擬摩擦接觸有效性，并揭示了其中涉及的許多細(xì)節(jié)。與基于庫侖摩擦和理想光滑接觸幾何形狀的傳統(tǒng)輪軌接觸模型相比，本模型描述了微米尺寸的粗糙表面的摩擦相互作用。更具體地說，該模型通過模擬軌道和車輪之間每個微凸體接觸對的生命周期來模擬摩擦接觸：從建立接觸到拉伸，完全接觸，到最后斷開并在短暫的休息時間后重新建立。它不僅可以描述輪軌接觸的縱向蠕變效應(yīng)，還可以揭示滑動速度對輪軌接觸摩擦性能的影響。

該研究所提出的模型，在滾動和滑動共同作用下，深刻地表征了輪軌接觸行為并揭示了摩擦演變的機制，這對于理解輪軌的動態(tài)接觸過程十分重要。

2 輪軌黏著力模型

采用的二維輪軌黏著力動態(tài)計算模型是基于Matlab程序設(shè)計的數(shù)值計算模型。考慮到粗糙度對輪軌接觸的影響，模型中利用微凸體之間的接觸來模擬輪軌間的蠕滑力。由于僅考慮縱向蠕滑，根據(jù)對稱性模型簡化為單輪單軌之間的接觸。

2.1 粗糙度的模擬

微觀尺度上，輪軌表面并非絕對光滑，而是有很多被稱為微凸體的不規(guī)則粗糙峰存在。兩個粗糙物體的接觸實質(zhì)上是微凸體之間的接觸，會導(dǎo)致真實接觸面積遠(yuǎn)小于名義接觸面積。在力學(xué)中，往往將兩個粗糙表面之間的接觸簡化為光滑表面和粗糙表面的接觸，故本文將輪軌接觸模擬為光滑車輪和粗糙鋼軌的接觸。結(jié)合Daves[15]等人實測的輪軌表面粗糙度的數(shù)據(jù)，假設(shè)輪軌表面的微凸體高度值的分布服從Gaussian分布，可以隨機生成一系列的微凸體數(shù)據(jù)，如圖1所示。

圖1 輪軌表面接觸示意

在數(shù)值模擬的過程中，在橢圓形的接觸內(nèi)隨機產(chǎn)生各個大小不等互不相交的圓來模擬微凸體之間的接觸。在本模型中，利用Matlab程序產(chǎn)生系列的微凸體，圖2即為數(shù)值模擬示意。圖中每個小圓代表一對相互接觸的微凸體接觸對，兩兩之間互不影響。因此在粗糙度的模擬中，微凸體的高度hi與半徑ri的關(guān)系為正比，即hi∝ri。兩個物體接觸時的輪廓接觸面積占名義接觸面積的5%～15%，根據(jù)這個比值可以大致確定接觸斑內(nèi)需要生成的微凸體接觸對的數(shù)量。圖2是通過Matlab程序在橢圓接觸內(nèi)隨機產(chǎn)生微凸體接觸圓的數(shù)值模擬結(jié)果。圖中的橢圓為接觸斑的輪廓，橢圓內(nèi)的小圓則代表著微凸體接觸對。

圖2 接觸斑微凸體接觸的數(shù)值模擬

2.2 接觸模擬

2.2.1 垂向接觸

將輪軌表面視作光滑的橢球面，利用Hertz接觸理論計算出接觸斑的幾何形狀和接觸內(nèi)的垂向應(yīng)力分布，如圖3所示，具體公式如下

(1)

(2)

(3)

式中，a、b分別為橢圓接觸斑的長短軸長度；P為單輪的軸重；pz(x,y)為接觸斑內(nèi)法向接觸應(yīng)力的分布；Rx1,Rx2,Ry1,Ry2分別為車輪和鋼軌接觸點處橫縱向曲率半徑；E(e)，e為橢圓積分參數(shù)。

圖3 輪軌垂向接觸示意

由于Hertz接觸計算的是無摩擦的光滑表面，而數(shù)值模擬中需要考慮到粗糙度，因此需要對法向應(yīng)力的數(shù)值進(jìn)行一定的修正，修正公式可表述為

H×(∑?i×si)=P

(4)

(5)

式中，αi為微凸體的受壓程度；si為微凸體完全接觸時受壓面積(即小圓面積)；H為材料硬度；P為單輪軸重。

σi為小圓中心處的壓應(yīng)力，由于ri?a(b)，可用小圓中心處的垂向應(yīng)力表示整個小圓上的平均應(yīng)力。

在不存在第三介質(zhì)(干摩擦)的情況下，輪軌相互作用的力全部由微凸體之間的接觸來傳遞。輪軌接觸時，并不是所有的微凸體之間都能達(dá)到預(yù)設(shè)的最大接觸面積，小圓面積si乘上受壓系數(shù)αi即表示當(dāng)前時刻的受壓面積。

2.2.2 切向接觸

切向接觸的模擬是本模型的重點，著重于從物理本構(gòu)的角度來描述摩擦。因此暫時只考慮縱向蠕滑率，只考慮單輪在單軌上的運動，簡化為橢圓接觸斑在生成的粗糙表面上移動。接觸斑的移動速度與車輪的行進(jìn)速度v0相同，車輪的微滑速度v1則使得彈簧拉伸來模擬蠕滑力。在微凸體相互接觸時，會經(jīng)歷兩個微凸體從相互接近到接觸，接觸面積逐漸到最大值之后接觸對逐漸分離的過程。半徑越大的微凸體在接觸時越難產(chǎn)生相對位移，即用來模擬接觸的彈簧上升剛度k1越大，即k1∝hi。在微凸體完全接觸(即彈簧到達(dá)受力極限)后，會在此狀態(tài)維持一段時間再開始分離，Ramin Aghababaei[14]利用離散元的方法模擬了整個微凸體從接觸到分離的全過程，為了簡化模型用線性彈簧來模擬整個受力過程。每個彈簧受力極限的確定公式可表述為

fi=μ·αi·si·H

(6)

式中，μ為垂向力和切向力的相關(guān)系數(shù)，在低速(微滑)的情況下，摩擦力仍與壓力呈正比。文獻(xiàn)[16]通過一系列的實驗確定了微凸體處于邊界潤滑時的摩擦系數(shù)，一般在0.2～0.4。當(dāng)模擬工況為干摩擦?xí)r，的取值可參考文獻(xiàn)[17]，一般為0.53～0.63。

圖4為單個彈簧的受力，兩個微凸體從開始接觸到完全接觸這一個階段，即彈簧拉力上升的階段。隨著滑動量的增大，微凸體接觸對之間的接觸面積隨之變大，切向力也隨著線性增大，直至兩個微凸體完全接觸。之后隨著車輪繼續(xù)滑動，微凸體接觸對開始分離，這時最大切向力會有一個短暫的穩(wěn)定期。隨著微凸體接觸對的分離，接觸力逐漸變小，直至完全分離?；赗amin Aghababaei的仿真數(shù)據(jù)和輪軌表面微凸體的幾何參數(shù)，微凸體分離時力的減小，假設(shè)彈簧下降段的剛度k2與k1正相關(guān)，k2∝k1。完全分離后的微凸體再次接觸之前會經(jīng)歷一個時間來恢復(fù)形變。從物理學(xué)的角度來看，兩個微凸體從相互分離到重新接觸，必然會經(jīng)過一個時間Δt。O.M.Barun的摩擦模型則驗證了Δt對描述摩擦的重要性以及物理意義，尤其是從靜摩擦到動摩擦的變化階段，這與接觸斑內(nèi)黏著區(qū)到滑動區(qū)的過度有著極高的相似性。

圖4 彈簧受力示意

當(dāng)前時刻的縱向蠕滑力即為接觸斑內(nèi)所有彈簧力的總和。在不同的時間步長內(nèi)，每個微凸體接觸對在接觸斑內(nèi)的位置發(fā)生改變，隨之而改變的是微凸體的法向受力，進(jìn)一步影響彈簧的受力極限fi，而不是嚴(yán)格如圖4所示，峰值會維持在一個恒定值。

2.3 模型概述

數(shù)值模擬利用MATLAB程序設(shè)計計算，計算流程如圖5所示。模型首先輸入車輪和鋼軌的參數(shù)，利用Hertz接觸理論算出單輪單軌工況下接觸斑的幾何參數(shù)和垂向應(yīng)力分布，經(jīng)過粗糙度修正之后利用彈簧摩擦模型進(jìn)行切向接觸模擬。

圖5 計算流程

3 結(jié)果分析

本文將通過模型所計算的牽引系數(shù)與蠕滑率的關(guān)系與試驗和線路測試的結(jié)果相比較，來驗證模型的正確性。本節(jié)將討論列車速度、縱向蠕滑率和粗糙度參數(shù)等因素對黏著系數(shù)的影響以及接觸斑內(nèi)滑動區(qū)和黏著區(qū)的分布。

3.1 模型的驗證

模型現(xiàn)階段主要是對干摩擦的模擬，考慮到對滾機難以模擬高速下的輪軌，現(xiàn)場所測得的高速列車數(shù)據(jù)也較少，故模型驗證選取中低速工況。圖6分別展示了列車速度為36 km/h[18]和60 km/h[19]時，仿真結(jié)果與現(xiàn)場測試以及對滾機實驗結(jié)果的對比，工況均為干摩擦。從圖6可以看出，數(shù)值模擬的結(jié)果與實測數(shù)據(jù)吻合程度很高。牽引系數(shù)與蠕滑率的曲線在小蠕滑時呈線性上升，達(dá)到峰值之后有一個下降的過程。

圖6 牽引系數(shù)-縱向蠕滑率曲線

表1顯示了模型在數(shù)值計算中使用的模擬參數(shù)。對于表格中尚未列出的參數(shù)如微凸體的最大靜摩擦系數(shù)則根據(jù)實際工況選取。

3.2 速度與蠕滑率

隨著高速鐵路的普及，而速度與蠕滑率對牽引系數(shù)影響非常大，因此研究牽引系數(shù)與之的關(guān)系十分重要。圖7展示了干摩擦工況下滾動速度、蠕滑率與牽引系數(shù)的關(guān)系。隨著速度的增加，黏著系數(shù)逐漸減小，這一結(jié)果符合現(xiàn)階段的認(rèn)識。隨著滾動速度的增大(接觸斑移動速度的增加)，在相同的時間內(nèi)，微凸體之間的相對位移更小，也即彈簧的拉伸量變小，蠕滑力因此變小。

表1 計算參數(shù)

圖7 牽引系數(shù)-速度-蠕滑率曲線

而在速度一定時，隨著蠕滑率的增大，牽引系數(shù)先增加后減小，與Kalker的接觸理論的計算結(jié)果相比，多出了之后的下降段。主要原因是后者在計算蠕滑力時，摩擦系數(shù)是一個常數(shù)，未考慮到動摩擦系數(shù)會隨著速度的增加而減小的現(xiàn)象。在速度一定時，隨著蠕滑率的變大，在相同的時間內(nèi)，彈簧的拉伸量變大。在彈簧的受力達(dá)到極限值之前，隨著拉伸量的增大，彈簧受力也變大，反映出牽引系數(shù)隨蠕滑率的變大而變大的規(guī)律，直至接觸斑到達(dá)完全滑動的狀態(tài)。當(dāng)縱向蠕滑率繼續(xù)增大時，切向接觸彈簧將會逐漸拉伸至彈簧的破壞極限值，即每個微凸體的接觸對開始分離，彈簧所表示的切向力開始下降，揭示了牽引系數(shù)-縱向蠕滑率曲線中下降段的原因。

3.3 粗糙度參數(shù)

輪軌表面由于加工等因素會不可避免地產(chǎn)生凹凸不平的微凸體。一般而言，表面粗糙度的關(guān)鍵參數(shù)主要是微凸體高度的均值和標(biāo)準(zhǔn)差即分布。本小節(jié)主要討論不同微凸體高度均值和標(biāo)準(zhǔn)差與均值之比對輪軌黏著特性的影響。

圖8為干摩擦工況下牽引系數(shù)隨著微凸體高度均值變化的曲線。從圖8可以看出，黏著系數(shù)是隨微凸體均值的增加而變大的。這種趨勢解釋了通過合理撒砂提高輪軌黏著系數(shù)的現(xiàn)象，結(jié)論和日本忠夫[20]試驗所測得的變化趨勢不一致，可能是由于濕摩擦?xí)r，隨著微凸體高度均值變大，接觸的微凸體總面積變小，而介質(zhì)承擔(dān)更多的接觸力導(dǎo)致黏著系數(shù)變??；干摩擦下仍由微凸體承擔(dān)全部的接觸力故未有變小的趨勢。

圖8 微凸起高度值均值-黏著系數(shù)曲線

由圖9可以看出，牽引系數(shù)隨著粗糙度分布的高度均值與標(biāo)準(zhǔn)偏差之比的增大而變大。在實際情況中，h/σ不會無限度的變大，具體取值可根據(jù)實地測量的輪軌粗糙度數(shù)值來確定。

圖9 標(biāo)準(zhǔn)差與高度均值之比-牽引系數(shù)曲線

輪軌間的黏著系數(shù)隨著表面粗糙度高度均值變大的主要原因是：隨著微凸體高度均值的增大，接觸區(qū)的相互接觸的微凸體由更多更小向更少更大轉(zhuǎn)變，即微凸體接觸對完全接觸時的接觸力更大，這可能是導(dǎo)致微凸體高度均值增大時，牽引系數(shù)也隨之增大的主要原因。從圖8也可以看出，在同一粗糙度下，車輪的滾動速度越大黏著系數(shù)越小，這種現(xiàn)象與3.2節(jié)中的結(jié)論一致。由此可見，某種程度上粗糙度是增加輪軌黏著系數(shù)的有利因素。在車輪踏面和鋼軌輪廓的設(shè)計和維護(hù)中，如果適當(dāng)?shù)目刂票砻娲植诙?，有利于改善黏著能力?/p>

3.4 摩擦系數(shù)

與傳統(tǒng)的輪軌力計算模型不同的是，模型中并未定義動摩擦系數(shù)，也未在計算中直接在接觸斑內(nèi)劃分滑動區(qū)和黏著區(qū)。模型通過計算接觸斑內(nèi)每個單元的牽引系數(shù)來劃分滑動區(qū)和黏著區(qū)。如圖10所示，在滾動速度為360 km/h時，隨著縱向蠕滑率的增大，接觸斑內(nèi)各個單元的牽引系數(shù)變化趨勢較為明顯。黃色的部分代表當(dāng)前單元的牽引系數(shù)近乎于最大值，即該單元內(nèi)的切向彈簧受力近似于最大值。在縱向蠕滑率小于0.6%時，達(dá)到最大牽引系數(shù)的單元由接觸斑末端向中心處逐漸增多。若劃定接觸斑內(nèi)的滑動區(qū)與黏著區(qū)，切向彈簧受力達(dá)到極限則微凸體的接觸對開始相對滑動，故可將牽引系數(shù)近似最大值的單元當(dāng)作滑動區(qū)。當(dāng)縱向蠕滑率繼續(xù)增大時，黃色單元會繼續(xù)向接觸斑前緣移動，滑動區(qū)的范圍進(jìn)一步擴大。同時由于接觸斑末端的接觸對開始相對滑動，切向彈簧受力值降低，導(dǎo)致接觸斑末端的單元牽引系數(shù)降低，但這些單元仍處在滑動區(qū)。

圖10 不同蠕滑率下的牽引系數(shù)分布

圖11是列車速度為360 km/h，縱向蠕滑率為0.8%時的牽引系數(shù)曲線。從圖11可以看出，牽引系數(shù)在接觸區(qū)前緣有一個線性增長的過程，對應(yīng)黏著區(qū)的摩擦系數(shù)是逐漸變大的，達(dá)到峰值時即為統(tǒng)計學(xué)上的最大靜摩擦系數(shù)，進(jìn)入到滑動區(qū)。牽引系數(shù)的驟降段是從最大靜摩擦系數(shù)到動摩擦系數(shù)的減小過程。在滑動區(qū)牽引系數(shù)繼續(xù)下降是由于微凸體接觸對開始相對滑動，切向接觸力逐漸減小。

圖11 牽引系數(shù)曲線(v=360 km/h，縱向蠕滑率ξ=0.8%，x軸正方向為列車前進(jìn)方向)

將牽引系數(shù)的峰值定義為本模型中的摩擦系數(shù)，可得切向力與垂向力乘摩擦系數(shù)的對比圖，如圖12所示。與傳統(tǒng)的蠕滑模型所不同的是，本模型采用的摩擦系數(shù)即非動摩擦系數(shù)，也非靜摩擦系數(shù)，而是通過統(tǒng)計學(xué)得到的一個平均靜摩擦系數(shù)，其值小于最大靜摩擦系數(shù)，大于動摩擦系數(shù)。由圖12的對稱性來看，此工況下接觸斑內(nèi)應(yīng)處于全滑動狀態(tài)，與圖10中縱向蠕滑率為0.8%的結(jié)果對應(yīng)，圖中藍(lán)色部分表示該單元內(nèi)此時刻沒有接觸體處于接觸狀態(tài)，接觸斑前緣有少量單元的牽引系數(shù)并未達(dá)到最大值，從統(tǒng)計學(xué)的角度來看該單元內(nèi)的大部分微凸體接觸對已經(jīng)滑動至接觸分離的狀態(tài)。

圖12 接觸力對比曲線(x軸正方向為列車前進(jìn)方向)

4 結(jié)論與展望

本文介紹一種中尺度動態(tài)輪/軌接觸模型，通過該模型，研究了影響牽引系數(shù)的輪軌速度，蠕變速率和粗糙度等因素。在分析數(shù)值計算結(jié)果之后，得到以下結(jié)論。

(1)當(dāng)考慮到粗糙度的影響時，接觸應(yīng)力的最大值比光滑表面的計算結(jié)果要大，且曲面不平滑存在著許多小鋸齒。

(2)隨著輪軌表面的微凸體高度值的增大，高速情況下的輪軌黏著系數(shù)變大。

(3)模型的計算結(jié)果與試驗數(shù)據(jù)吻合良好，驗證了數(shù)值計算模型和算法的可行性。這表明該模型雖然被設(shè)計為動態(tài)模型，但仍能夠預(yù)測穩(wěn)態(tài)牽引因子。

(4)模型中并未人為地劃定滑動區(qū)與黏著區(qū)，而是通過整個接觸斑內(nèi)各個單元牽引系數(shù)的變化，刻畫出滑動區(qū)到黏著區(qū)過渡的整個變化過程。

在以后的研究中會將模型從二維拓展到三維，進(jìn)一步考慮到橫向蠕滑和自旋蠕滑的影響，并進(jìn)行高速下的模型驗證。