彈性半空間非赫茲輪軌法向壓力分布

車輪接觸滑動度的計算是風-列車-橋梁耦合振動分析的中心問題之一。目前在鐵路部門應用較為普遍的輪軌接觸蠕滑力計算模型均為赫茲接觸型,如Kalker線性接觸理論、Kalker簡化理論和Shen-Hedrick-Elkins理論等。赫茲法向接觸理論認為接觸點處的曲率半徑為常數,接觸斑為橢圓,接觸斑內的壓力為半橢球體分布。實際上輪軌滾動接觸問題不完全符合赫茲接觸條件,如輪軌接觸型面是由多條不同曲率半徑圓弧組成,輪軌滾動接觸過程的接觸點常常位于這些不同半徑圓弧的交界處,這就不滿足接觸點處曲率半徑是常數的條件,輪軌之間彈性擠壓量和法向壓力分布的大小和接觸斑區域與按照赫茲理論求得的結果均存在一定的差異,在輪緣發生單接觸斑內兩點接觸時的這一差異尤其顯著。在列車-橋梁耦合振動分析中引入橫向風作用后會使得車輛相對于軌道發生較為明顯的橫移,從而使得輪軌發生輪緣接觸的頻率大為增加。在輪緣接觸情況下,赫茲接觸橢圓面積小于非赫茲接觸面積,最大赫茲正壓力大于最大非赫茲正壓力,這種差異是顯著的,此時輪軌之間還存在大自旋蠕滑,而目前常用的輪軌蠕滑力模型在計算大自旋蠕滑時均存在較大不足,基于赫茲接觸理論的蠕滑力計算理論僅能在工程精度內處理非輪緣單點接觸的情況。考慮到基于Kalker精確理論的DUVOROL數表和TPLP數表并未公開發布,本文在以非赫茲接觸理論計算輪軌接觸法向壓力分布的基礎上,根據修正的FastSim程序計算了輪軌在單點接觸、輪緣接觸和單接觸斑內兩點接觸情況下的蠕滑力,并編制了便于風-列車-橋梁耦合分析應用的蠕滑力插值數表MFTTLM(ModifiedFastSimTractionTable-LM)。1接觸區域的離散化輪軌接觸非赫茲法向壓力分布的計算以彈性半空間理論為基礎,使用數值計算方法得到壓力的分布。兩個彈性體發生法向接觸時的控制方程為{w1+w2+f12=δ(x,y)∈Cw1+w2+f12＞δ(x,y)?Cp(x,y)≥0(x,y)∈Cp(x,y)=0(x,y)?C∫∫Cp(x,y)dxdy=Ρ(1)???????????????????????????w1+w2+f12=δw1+w2+f12＞δp(x,y)≥0p(x,y)=0∫∫Cp(x,y)dxdy=P(x,y)∈C(x,y)?C(x,y)∈C(x,y)?C(1)式中,C為接觸區域;w1、w2分別表示接觸切向平面中坐標同為(x,y),位于兩個彈性體1、2上的兩個點發生的彈性變形;f12為接觸切向平面(x,y)坐標處的法向間隙;δ為兩個彈性體沿接觸區域的法向分別在兩個無窮遠處的相對接近量;p(x,y)為法向壓力分布;P為總的法向壓力大小。對接觸區域C進行離散化,令每個單元內各點的彈性位移、法向間隙和壓力大小均為常數并以單元形心處的值wi、fi和pi來表達。考慮到車輪與鋼軌為兩個具有相同彈性常數(彈性模量E和泊松比ν)的接觸體,在每一個單元形心,控制方程(1)中第1式變形為wi=12(δ-fi)i=1,2,?,n(2)wi=12(δ?fi)i=1,2,?,n(2)根據Boussinesq關于彈性半空間在表面集中力作用下的應力與變形關系,在每一個單元形心均存在wi=n∑j=1φijpji=1,2,?,n(3)wi=∑j=1nφijpji=1,2,?,n(3)因此得到每個單元的位移方程n∑j=1φijpj=12(δ-fi)i=1,2,?,n(4)∑j=1nφijpj=12(δ?fi)i=1,2,?,n(4)式中,φij為影響系數,表示j單元內部的單位法向均布壓力在i單元形心處產生的法向位移。對于矩形單元,其表達式為φij=1-ν2πExj+aj∫xj-ajyj+bj∫yj-bjdxdy√(x-xi)2+(y-yi)2(5)對于2a×2b的矩形單元,Love給出了影響系數的解析表達式。應當注意到,在計算i=j的影響系數時,式(5)和Love給出的表達式均為奇異。考慮到在對接觸區域進行離散時,矩形單元的長寬比越大越容易產生壓力分布的奇異性,正方形單元在矩形類型單元中可以保持最好的離散精度,本文采用正方形對接觸區域進行離散,此時i=j情況下的影響系數可以通過對式(5)進行積分變換得到φij=1-ν2πE8aln(√2+1)(6)根據每個單元上的位移方程式(4)可以得到離散后的整個接觸區域的位移方程組[φ11φ12?φ1nφ21φ22?φ2n???φn1φn2?φnn〗[p1p2?pn〗=[(δ-f1)/2(δ-f2)/2?(δ-fn)/2〗記為[Φ][Ρ]=[D](7)位移方程式(7)中方程左邊的法向壓力pi與方程右邊的接近量δ均為未知量,可以使用迭代的方法進行求解,具體求解步驟為:①對可能的接觸區域進行離散,根據輪軌接觸幾何參數和輪軌踏面的外型計算接觸區域各個單元形心處的法向間隙fi(i=1,2,3,…,n)。②指定初始接近量δ,計算可能的接觸區域C。③形成當前接觸區域C的影響系數矩陣,求解法向壓力大小pi(i=1,2,3,…,n)。④判斷接觸區域C是否全部滿足控制方程(1)中第3式。如果否,則修正接觸區域C,進入步驟③。⑤判斷接觸區域C是否在一定的收斂誤差范圍內滿足控制方程(1)中第5式。如果否,則修正接近量δ,進入步驟②,否則計算結束。2輪軌接觸的彈性常數在計算得到輪軌接觸的接近量和法向壓力分布后,使用FastSim算法計算輪軌接觸區域上的蠕滑力。FastSim算法以Kalker的簡化理論為基礎,但并不局限于橢圓接觸斑,還可以應用于多點接觸的情況并可以得到較好的結果。當將FastSim算法應用于非橢圓接觸區域中時,關鍵的問題在于合適地確定接觸區域的彈性常數L。彈性常數L與接觸等效橢圓的半軸長ai、bi和半軸長之比ai/bi相關。當輪軌發生單接觸斑內單點接觸時,將非赫茲法向接觸計算得到的接觸區域使用長度為l、寬度為w的最小矩形面積包容,定義λ=lw(8)作為對赫茲接觸橢圓半軸長比的修正,同時根據下式對赫茲接觸橢圓進行尺寸修正{aebe=λπaebe=A(9)式中,ae、be為等效接觸橢圓的尺寸;A為非赫茲接觸區域的面積。以此便可以計算彈性常數L。在得到相對赫茲接觸情況進行修正之后的彈性常數,在真實的接觸區域之上根據真實的法向壓力分布應用FastSim算法,便可計算得到輪軌接觸的縱向和橫向蠕滑力。對于輪軌接觸單接觸區域內兩點接觸的情況,利用上面的彈性常數確定方法對不同的等效橢圓確定不同的彈性常數,在每一個等效接觸橢圓所對應的真實接觸區域上應用FastSim算法,計算得到輪軌接觸的縱向和橫向蠕滑力。根據上述的輪軌蠕滑力計算方法(修正的FastSim),考慮LM磨耗型踏面和60kg/m鋼軌,軌底坡為1∶40,輪緣內側距為1353mm,軌距為1435mm,名義滾動圓半徑r0=420mm,輪重P0=105kN,鋼軌和車輪為同質材料,取彈性模量E=210GPa,泊松比ν=0.28,靜摩擦系數fs=0.3,動摩擦系數fk=0.29,取行車速度v=30m/s,搖頭角位移ψ和速度˙ψ均為0,輪對橫移y=0~12mm?˙y/ν=0、0.0001、0.0005和0.0010分別對應工況(1)、(2)、(3)和(4),計算得到的左右輪對接觸點處的縱向和橫向蠕滑力隨輪對橫移的變化見圖1和圖2。3輪緣接觸的模型模擬Kalker的三維彈性體非Hertz滾動接觸理論(精確理論)是到目前為止最完善的三維滾動接觸理論,基于此理論的計算機程序CONTACT可以求解非Hertz接觸條件下的輪軌滾動接觸蠕滑問題。考慮到修正的FastSim著重于處理輪緣接觸時的蠕滑力求解,以下僅列出修正的FastSim與CONTACT在計算輪緣接觸時兩種蠕滑力結果的對比見表1。計算輪軌蠕滑力時采用TB錐形踏面和60kg/m鋼軌,軌底坡為1∶40,輪緣內側距為1353mm,軌距為1435mm,名義滾動圓半徑r0=420mm。接觸斑內總的法向壓力N=83.6kN,鋼軌和車輪為同質材料,取剪切模量G=82GPa,泊松比ν=0.28,靜摩擦系數fs=0.3,動摩擦系數fk=0.29。輪對的橫移為11.0mm,搖頭角位移為1.0°。計算所采用的蠕滑率數據確保了輪緣接觸時接觸斑經歷由負向滑動(飽和)→完全黏著→正向滑動(飽和)的整個過程,具有較強的代表性。由表1可見,修正的FastSim與CONTACT的縱向、橫向蠕滑力計算誤差不超過5%,蠕滑力的合力誤差小于1%,具有比較精確的計算結果。因此可以認為,修正的FastSim在計算輪緣接觸時的蠕滑力是可靠的。LM磨耗型踏面和60kg/m鋼軌在輪對橫移為10.0mm、搖頭角位移為0°時發生單接觸斑內兩點接觸,表2為不同蠕滑率下的各種蠕滑力計算方法得到的蠕滑力差異,計算條件如本節開始時所述。由表2可知,在計算兩點接觸時,Shen-Hedrick-Elkins(S-H-E)理論計算得到的蠕滑力與修正的FastSim(MF)計算得到的蠕滑力相比較,縱向蠕滑力的最大相對誤差為11.06%,橫向蠕滑力的最大相對誤差為-34.48%。而FastSim與修正的FastSim計算結果對比表明,縱向蠕滑力的最大相對誤差為-7.73%,橫向蠕滑力的最大相對誤差為11.87%。考慮到對于輪軌接觸單接觸區域兩點(或多點)接觸的情況,分別考慮各個接觸等效橢圓計算蠕滑力的效果要好于將輪軌接觸視為單點接觸,可以認為修正的FastSim算法在計算單接觸斑內兩點接觸的蠕滑力時的精度相對于Shen-Hedrick-Elkins理論和FastSim算法均有較大的提高。4平均速度計算盡管修正的FastSim在計算輪緣接觸和單接觸斑內兩點接觸情況下的蠕滑力時相對目前通用于列車-橋梁耦合中的蠕滑力計算方法可以給出更加準確的結果,更加適合應用于風-列車-橋梁耦合振動分析之中以便更好地反映橫向風對列車-橋梁耦合的影響,但是對于風-列車-橋梁耦合動力時程分析所遭遇的大規模數值迭代分析而言,修正的FastSim算法計算蠕滑力的速度(在Pentim4,2.66GHz的CPU,512M的DDR內存配置下,計算一次蠕滑力的平均時間為6s)是不可接受的。因此,針對當前鐵路車輛車輪踏面的行業標準,基于修正的FastSim編制了LM磨耗型踏面與60kg/m鋼軌匹配時的蠕滑力插值數表MFTTLM(ModifiedFastSimTractionTable-LM)。編制插值數表的基本原理見文獻。5接觸模式分析考慮橫向風影響之后的列車-橋梁耦合振動分析需要準確地處理更多幾率的輪緣接觸,同時為了更好地處理單接觸斑內兩點接觸的情況,在以非赫茲接觸理論計算輪軌接觸法向壓力分布的基礎上,根據修正的FastSim算法計算