重載鐵路鋼軌磨耗狀態(tài)下的輪軌法向接觸特性

許玉德，嚴(yán)道斌，孫小輝，唐永康

(1.同濟(jì)大學(xué) 道路與交通工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 201804；2.上海市軌道交通結(jié)構(gòu)耐久與系統(tǒng)安全重點(diǎn)實(shí)驗(yàn)室，上海 201804；3.比亞迪汽車工業(yè)有限公司，廣東 深圳 518118；4.朔黃鐵路發(fā)展有限責(zé)任公司，河北 滄州 062356)

重載列車速度和軸重的不斷提高，加劇了重載鐵路鋼軌的磨耗.磨耗使得輪軌接觸狀態(tài)發(fā)生改變，輪軌接觸由單點(diǎn)接觸變?yōu)閮牲c(diǎn)接觸或共形接觸的概率大大提升，輪軌接觸界面的作用變得更加復(fù)雜.研究表明，鋼軌磨耗與輪軌法向接觸應(yīng)力直接相關(guān)[1-2].

國內(nèi)外學(xué)者對輪軌接觸理論展開了深入研究.Kalker[3-5]提出三維彈性體滾動接觸理論，基于此開發(fā)的CONTACT程序，可以求解考慮蠕滑的輪軌接觸問題；為了適用車輛動力學(xué)仿真，Kalker提出簡化理論并開發(fā)了FASTSIM程序.Li[6-7]對Kalker在CONTACT中采用的彈性半空間影響系數(shù)作了修正，提出四分之一空間的影響系數(shù)有限元計算方法，同時在計算輪軌接觸時不再將蠕滑視作常數(shù)，基于此開發(fā)的WEAR程序，能夠有效處理軌距角輪軌共形接觸問題.金學(xué)松[8]在Kalker理論的基礎(chǔ)上，分析了輪軌滾動接觸蠕滑率與蠕滑力之間的關(guān)系，拓展了其簡化理論模型，使其更適合求解輪軌黏著問題.任尊松[9]在跡線法的基礎(chǔ)上提出了跡線極值法，建立車輛-軌道系統(tǒng)動力學(xué)模型，能夠解決輪軌多點(diǎn)接觸問題并獲得準(zhǔn)確的輪軌多點(diǎn)接觸幾何參數(shù).楊新文[10]提出基于赫茲理論的P_M(partition-model)模型，采用切片投影法尋求接觸區(qū)域，并利用分區(qū)方法計算非橢圓接觸斑上的法向接觸應(yīng)力.

當(dāng)鋼軌產(chǎn)生磨耗時，車輪與鋼軌的接觸可能會發(fā)生在軌距角區(qū)域，尤其是在小半徑曲線上.軌距角區(qū)域的輪軌接觸是曲面接觸，彈性半空間假設(shè)不再成立[11].CONTACT程序由于采用彈性半空間假設(shè)，因此不能較好地處理這類問題；而有限元方法不受彈性半空間假設(shè)的限制，計算結(jié)果準(zhǔn)確，是目前處理軌距角接觸問題常用的選擇[12]，其缺點(diǎn)是計算效率較低.本文基于Kalker三維彈性體滾動接觸理論，引入法向角對彈性半空間假設(shè)影響系數(shù)進(jìn)行修正，在局部坐標(biāo)系下計算輪軌接觸的法向間隙，對修正后的Kalker法向接觸方程編程求解，可以準(zhǔn)確、快速地計算輪軌法向接觸應(yīng)力，提出的方法適用于輪軌平面和非平面接觸.

1 輪軌法向接觸應(yīng)力計算

1.1 輪軌滾動接觸理論

Kalker的三維彈性體滾動接觸理論，是一種基于邊界元的數(shù)值計算方法.如圖1所示，該方法將可能的接觸區(qū)域劃分為矩形單元網(wǎng)格，利用彈性半空間Boussinesq-Cerruti公式計算得到任一單元與可能接觸區(qū)域內(nèi)其他所有單元的位移差，借助Newton-Raphson方法對接觸方程迭代求解得到所有接觸單元的應(yīng)力分布[13].

圖1 可能接觸區(qū)域的網(wǎng)格劃分Fig.1 Meshing of possible contact region

對于法向接觸問題，Kalker理論中最小余能方程為

(1)

式中：C為余能；I、J為矩形單元的編號；pIz為作用在I單元的法向應(yīng)力；AIzJz為影響系數(shù)，表示作用在J單元的法向單位力在I單元處引起的法向位移；hI為I單元處輪軌的法向間隙；δ為法向侵入量；S0為每個矩形單元的面積，S0=Δx1·Δx2；M為接觸區(qū)域內(nèi)單元數(shù)量；Fz為法向力.

1.2 影響系數(shù)的修正

影響系數(shù)AIzJz對法向接觸方程的求解精度起到重要作用，在法向接觸問題研究中影響系數(shù)的計算多采用Boussinesq-Cerruti彈性半空間解析解.當(dāng)鋼軌產(chǎn)生磨耗時，車輪與鋼軌可能在軌距角區(qū)域接觸，彈性半空間假設(shè)不再成立.有限元計算影響系數(shù)的方法不受彈性半空間假設(shè)限制，因此計算影響系數(shù)的結(jié)果是準(zhǔn)確的.文獻(xiàn)[14]利用有限元方法計算影響系數(shù)，但計算效率較低.文獻(xiàn)[6]認(rèn)為彈性半空間假設(shè)在處理法向接觸問題時仍然是一種合理的選擇.文獻(xiàn)[15]提出了一種考慮法向角的非平面影響系數(shù)計算方法.

如圖2所示，將鋼軌視作二維平面內(nèi)的曲線，不考慮縱向方向的影響，計算J單元受法向單位力作用下可能接觸區(qū)域內(nèi)I單元的影響系數(shù)，處理過程遵循以下步驟：

(1)擬合鋼軌型面，選取可能接觸區(qū)域內(nèi)I單元，確定其法向方向nI和切向方向sI；

(2)對作用在J單元的法向單位力F按照nI和sI分解為FnI和FsI；

(3)設(shè)I單元和J單元的切線斜率分別為kI和kJ，以確定I單元和J單元的法向夾角α.

圖2 影響系數(shù)計算Fig.2 Calculation of influencing coefficient

完成(1)～(3)步驟后，可以近似認(rèn)為I單元在某一個平面內(nèi)受到一個法向力FnI和一個切向力FsI的作用，由于是在局部坐標(biāo)系中的計算，故局部坐標(biāo)系的影響系數(shù)可以表示為

AInJn=AIzJzcosα-AIzJysinα

(2)

式中：AInJn為局部坐標(biāo)系的影響系數(shù)；AIzJz和AIzJy為全局坐標(biāo)系的影響系數(shù)，由Boussinesq-Cerruti公式計算得到；α為法向角.

以接觸區(qū)域的曲率半徑15 mm為例，進(jìn)行影響系數(shù)計算.選取CHN75鋼軌型面，在軌距角處劃定一個10 mm×10 mm的矩形可能接觸區(qū)域，每個單元是0.5 mm×0.5 mm的正方形單元，在矩形區(qū)域的中心點(diǎn)處作用一個單位法向力，計算可能接觸區(qū)域內(nèi)各單元的影響系數(shù).為驗(yàn)證該近似方法的準(zhǔn)確性，利用有限元方法和CONTACT計算在相同條件下的影響系數(shù)，結(jié)果如圖3所示.

圖3 影響系數(shù)計算方法對比Fig.3 Comparison of calculation methods of influencing coefficient

從影響系數(shù)的最大值來看，有限元方法計算所得影響系數(shù)最大值為3.06×10-6mm，本文方法為3.10×10-6mm，誤差為+1.31%；CONTACT基于Boussinesq-Cerruti彈性半空間解析解的結(jié)果為2.92×10-6mm，誤差為-4.6%.從變化趨勢來看，3種方法計算的結(jié)果趨勢比較一致，但本文方法與有限元方法得到的結(jié)果更接近，因此本文采用的方法是準(zhǔn)確的.

1.3 法向間隙的修正

法向間隙是計算法向接觸問題的重要輸入?yún)?shù).常用的輪軌法向間隙計算方法是將輪軌接觸看作是全局坐標(biāo)下的接觸，然后計算兩點(diǎn)之間沿y軸方向的差值，作為法向間隙.CONTACT程序中認(rèn)為任一單元的法向間隙可用拋物線形函數(shù)表征

h(x，y)=b1x2+b2xy+b3y2+b4x+b5y+b6

(3)

式中：(x，y)為單元的全局坐標(biāo)；b1、b2、…、b6均為根據(jù)車輪和鋼軌型面確定的參數(shù).

這種方法仍然是基于彈性半空間假設(shè)，其認(rèn)為鋼軌接觸表面為平面或者近似平面，輪軌在軌距角接觸時顯然是不成立的.本文采用在局部坐標(biāo)系中的一種近似處理的方法，如圖4所示，將車輪與鋼軌視作二維平面內(nèi)的兩條曲線，不考慮縱向的影響，計算可能接觸區(qū)域內(nèi)P單元的法向間隙，處理過程遵循以下步驟：

(1)擬合車輪與鋼軌型面，對鋼軌表面進(jìn)行網(wǎng)格劃分，選取可能接觸區(qū)域內(nèi)的P單元，作P單元的法線，交車輪型面曲線于Q0點(diǎn)；

圖4 法向間隙的計算Fig.4 Calculation of normal distance

為驗(yàn)證這種方法的準(zhǔn)確性，采用1.2節(jié)中的計算參數(shù)，補(bǔ)充的參數(shù)為選取LMA系列原始車輪型面，輪對橫移量為+12 mm，θ取0.1°，Δ取10-5mm，利用有限元方法和CONTACT計算在相同條件下的法向間隙，結(jié)果如圖5所示.

圖5 法向間隙計算結(jié)果對比Fig.5 Comparison of normal distance results

本文方法與有限元方法的計算結(jié)果比較接近，尤其是在接觸點(diǎn)中心附近區(qū)域.從利用CONTACT拋物線多項式差值計算的結(jié)果來看，這種方法在處理非平面問題時存在較大誤差，尤其是在離接觸點(diǎn)中心較遠(yuǎn)的地方.在橫向方向-5 mm處有限元方法法向間隙的結(jié)果為0.032 mm，本文方法計算結(jié)果為0.031 mm，CONTACT計算結(jié)果則為0.051 mm，相比于有限元結(jié)果誤差達(dá)59.4%；在+5 mm處有限元方法法向間隙的結(jié)果為0.050 mm，本文方法計算結(jié)果為0.049 mm，CONTACT計算結(jié)果則為0.071 mm，相比于有限元結(jié)果誤差達(dá)42.0%.本文方法與有限元方法結(jié)果較為一致，因此準(zhǔn)確性較好.

基于上述影響系數(shù)和法向間隙的修正，對Kalker三維彈性體滾動接觸理論中法向接觸問題的最小余能方程進(jìn)行編程求解，可以準(zhǔn)確、快速地計算輪軌接觸的法向應(yīng)力分布情況.

2 鋼軌磨耗狀態(tài)下的法向接觸特性

以30 t軸重列車為例進(jìn)行計算，車輪選取LMA系列原始型面，鋼軌選取CHN75軌通過總重達(dá)100 Mt的磨耗型面，如圖6所示，該磨耗型面測量于我國某條重載鐵路.需要說明的是，本文后續(xù)圖中，均以鋼軌“工”字形截面的水平向?yàn)闄M向方向，鉛垂向?yàn)榇瓜蚍较颍劁撥夐L度方向?yàn)榭v向方向.

圖6 重載鐵路CHN75磨耗鋼軌Fig.6 Heavy haul CHN75 worn rail

2.1 軌頂面接觸

當(dāng)車輪的橫移量為0時，輪軌接觸為軌頂面接觸，屬于平面接觸范疇，計算結(jié)果如圖7所示.此時輪軌接觸斑為近似橢圓形，接觸面積為196.83 mm2，最大法向應(yīng)力為1 757 MPa.

2.2 軌距角接觸

當(dāng)輪對橫移量為+13 mm時，輪軌接觸為軌距角接觸，此時為兩點(diǎn)接觸，計算結(jié)果如圖8所示.此時接觸斑形狀不再是橢圓形，兩個接觸斑的面積分別為69.66 mm2和165.24 mm2，對應(yīng)的最大法向接觸應(yīng)力分別為1 582 MPa和1 740 MPa.

圖7 橫移量為0時的輪軌法向接觸應(yīng)力Fig.7 Normal contact stress at 0 mm lateral displacement

圖8 橫移量為+13 mm時的輪軌法向接觸應(yīng)力Fig.8 Normal contact stress at +13 mm lateral displacement

2.3 不同橫移量下的法向接觸特性

由于列車運(yùn)行時是蛇形運(yùn)動，因此輪對存在一定的橫移量，尤其在曲線運(yùn)行時，更是如此.利用跡線法確定在此鋼軌磨耗狀態(tài)下不同橫移量時的輪軌接觸點(diǎn)，選取-12～+14 mm的橫移量，每隔2 mm取一種工況，計算得到不同橫移量下的輪軌接觸應(yīng)力和接觸面積，見圖9.可以看出，在橫移量為-12～+10 mm時，接觸斑形態(tài)變化不大；在+12～+14 mm時，輪軌接觸由車輪踏面與鋼軌軌頂面接觸變?yōu)檐囕嗇喚壓弯撥壾壘嘟堑慕佑|，最大接觸應(yīng)力達(dá)到3 287 MPa.

圖9 不同橫移量下的法向接觸應(yīng)力Fig.9 Normal contact stress at each lateral displacement

最大法向應(yīng)力和接觸面積的計算結(jié)果如圖10所示.最大法向應(yīng)力的變化趨勢與接觸面積的變化趨勢相反，隨著接觸面積的減小，最大法向應(yīng)力相應(yīng)增大.

圖10 不同橫移量下最大法向應(yīng)力與接觸面積Fig.10 Maximum contact stress and contact area at each lateral displacement

在橫移量為+12.9～+13.2 mm區(qū)間時出現(xiàn)了兩點(diǎn)接觸的現(xiàn)象，如圖11所示.

圖11 不同橫移量下兩點(diǎn)接觸的法向接觸應(yīng)力Fig.11 Normal contact stress at two-point contact

對兩點(diǎn)接觸的兩個接觸斑進(jìn)行分析，最大法向應(yīng)力和接觸面積的結(jié)果如圖12和圖13所示.在兩點(diǎn)接觸的情況下，最大法向應(yīng)力與接觸面積變化的趨勢是一致的.

圖12 靠近鋼軌中心線的接觸斑最大法向應(yīng)力和接觸面積Fig.12 Maximum contact stress and contact area of the region closer to the rail center line

3 結(jié)論

基于Kalker法向接觸的最小余能方程，對影響系數(shù)和法向間隙進(jìn)行修正，引入法向角計算非平面接觸時的影響系數(shù)，在局部坐標(biāo)系下獲得接觸點(diǎn)的法向間隙，所采用的方法與有限元方法得到的結(jié)果接近，驗(yàn)證了方法的準(zhǔn)確性.對修正后的最小余能方程進(jìn)行編程，可以求解輪軌接觸的應(yīng)力分布，并選取某重載鐵路通過總重達(dá)100 Mt的CHN75型面磨耗鋼軌，車輪選取LMA系列原始型面，利用所編程序，計算了30 t軸重下重載磨耗鋼軌的輪軌法向接觸特性，得到以下結(jié)論：

圖13 靠近軌道中心線的接觸斑最大法向應(yīng)力和接觸面積Fig.13 Maximum contact stress and contact area of the region closer to the track center line

(1)當(dāng)橫移量為-12～+12 mm時，輪軌接觸為單點(diǎn)接觸，輪軌接觸狀態(tài)接近，最大法向應(yīng)力和接觸面積小幅變化，且最大法向應(yīng)力與接觸面積的變化趨勢相反.

(2)當(dāng)橫移量為+12.9～+13.2 mm時，輪軌接觸狀態(tài)為兩點(diǎn)接觸，在兩個接觸斑上最大法向應(yīng)力和接觸面積變化幅度較大，且最大法向應(yīng)力與接觸面積的變化趨勢相同.

(3)當(dāng)橫移量達(dá)到+14 mm時，輪軌接觸為車輪輪緣和鋼軌軌距角接觸，此時接觸斑變得狹長，接觸面積較小，但最大接觸應(yīng)力達(dá)到3 287 MPa.

本文對鋼軌磨耗狀態(tài)下的法向接觸特性進(jìn)行了研究，但輪軌法向接觸應(yīng)力與鋼軌磨耗之間的關(guān)系還有待更深入探討，這是筆者未來研究的重點(diǎn).同時，如何讓計算結(jié)果更適用于車輛動力學(xué)仿真也是需要進(jìn)一步研究的方向.