伏培林 丁立 趙吉中 張旭 闞前華,2) 王平

*(西南交通大學力學與工程學院,成都 610031)

?(西南交通大學土木工程學院,成都 610031)

引言

近年來,隨著我國軌道交通行業(yè)的迅速發(fā)展,列車的行駛速度和運量均不斷提升,這對列車的牽弓和制動性能提出了越來越高的要求.列車的牽弓和制動性能通過輪軌之間的接觸關系來實現(xiàn),一旦列車的牽弓力或制動力超過輪軌接觸關系所能提供的粘著力時,輪軌間便會發(fā)生相對滑動[1].在滑動過程中,輪軌接觸斑內的溫度會因摩擦生熱而迅速上升,并隨著車輪的進一步滑動脫離接觸而快速衰減,導致熱應力出現(xiàn)交變變化.這種交變的溫升誘發(fā)熱應力,疊加在軸重弓起的機械應力上,將大幅度降低鋼軌和車輪的服役壽命[2-3].同時,當滑動溫升足夠高時,輪軌表層材料發(fā)生馬氏體相變,并伴隨著后續(xù)的快速冷卻過程而形成脆硬的馬氏體組織.在輪軌接觸載荷的作用下,這些脆硬組織容易發(fā)生脫落,造成輪軌表面的扁疤和材料剝離,進而影響列車行駛的安全性[4-5].另外,過高的輪軌表面溫度也會導致材料硬度等力學性能的顯著下降,從而對輪軌表面的耐磨損性能帶來不利影響[6].由于對輪軌滑動接觸溫升進行實驗測量往往比較困難,因此通過數(shù)值模擬和解析分析等非測試方式開展相應的溫升研究具有重要的工程意義.

隨著ABAQUS、ANASYS 等商業(yè)有限元軟件的普及,輪軌滑動接觸溫升的有限元模擬研究受到了學者們的廣泛關注.然而,由于輪軌接觸斑非常狹小,為了準確地模擬接觸斑內的溫度場,在有限元模擬過程中必須對這一區(qū)域進行大規(guī)模網格加密處理,從而所花費的時間成本較為高昂.因此,針對輪軌滑動接觸溫升問題的解析或半解析研究也得到了大量開展.例如,Tanvir[7]忽略了傳熱參數(shù)的溫度相關性,將Hertz 接觸壓力的橢圓型分布函數(shù)近似成四階多項式級數(shù)形式,通過拉普拉斯變換給出了二維輪軌滑動摩擦溫升分布的解析解.基于相同的材料屬性溫度無關性假設,Knothe 和Liebelt[8]通過拉普拉斯變換和格林函數(shù)法,給出了任意接觸壓力下的輪軌滑動溫升問題的解析解,討論了法向接觸壓力分布形式和表面粗糙度對鋼軌表面溫度的影響.基于Hertz 彈性接觸理論,Ertz 和Knothe[9]討論了輪軌間換熱系數(shù)等因素對輪軌摩擦溫升的影響.Fischer 等[10]同時考慮了摩擦功和局部塑性變形對接觸界面熱源的貢獻,給出了輪軌滑動接觸閃溫的解析表達式.此外,考慮到輪軌間存在磨屑等雜物,Wu 等[11]基于三體微接觸模型和接觸溫升理論,討論了雜物尺寸、雜物密度、表面粗糙度、相對速度以及法向載荷等參數(shù)對接觸界面溫度場的影響.

由于材料屬性的溫度相關性,摩擦副材料的熱力學性質在摩擦溫升過程中也會發(fā)生變化,考慮材料屬性的溫度相關性將進一步提升溫升計算結果的準確性.然而,材料屬性溫度相關性的弓入會導致熱傳導方程的非線性化,大大增加了溫升求解的困難性.因此,考慮材料屬性的溫度相關性,必須首先求解非線性熱傳導方程.Ling 和Rice[12]考慮了熱導率和比熱容的溫度相關性,通過傅里葉余弦變換,給出了表面具有移動熱源半無限大體的準穩(wěn)態(tài)溫度場的迭代數(shù)值解法.這種求解方法允許任意形式的溫度相關性函數(shù),但相應的計算成本高,在實際應用中受到一定的限制.梁鈺等[13-14]基于特征正交分解方法,對非線性瞬態(tài)熱傳導問題的有限元離散格式進行降階處理,從而提出了一種求解考慮傳熱參數(shù)溫度相關性的二維非線性瞬態(tài)熱傳導問題的高效算法.通過假設組份材料和孔隙填充物的熱導率均隨溫度呈三次函數(shù)變化,戴婷等[15]采用微分求積法求解了一類含孔隙變厚度功能梯度圓盤的穩(wěn)態(tài)溫度場.當熱擴散系數(shù)隨溫度線性變化時,Fabre 和Hristov[16]基于熱平衡積分法給出了一類一維熱傳導方程的解析解.林府標等[17]通過李群分析法,考慮擴散率關于溫度的冪函數(shù)形式,得到了該類一維廣義熱傳導方程的顯示解析解和行波解.Ghasemi 等[18]則假定熱導率是溫度的線性函數(shù),將其斜率作為擾動小量,也給出了一類一維非線性熱傳導問題的攝動解法.此外,微分變換法[19]、有限差分法[20-21]和近場動力學法[22]也被應用于非線性熱傳導溫度場的求解.

由上述分析可見,已有研究存在不足:(1) 接觸壓力通常被假定為滿足Hertz 彈性接觸理論中的半橢圓型分布形式,而在實際輪軌接觸載荷作用下,接觸斑內將發(fā)生塑性變形,其接觸壓力不再滿足半橢圓型分布形式,必須進行塑性修正[23-25].(2) 材料參數(shù)往往被簡化為溫度無關的常量,或僅考慮某一個材料參數(shù)的溫度相關性.實際上,包含摩擦系數(shù)在內的多種材料屬性均會隨溫度的改變而發(fā)生變化[26],在溫升計算模型中同時考慮多種材料屬性的溫度相關性,有望提升溫升計算結果的準確性.

因此,本文通過擬合實驗數(shù)據(jù),獲得了熱導率、比熱容以及摩擦系數(shù)的溫度相關性函數(shù),并基于彈塑性接觸理論,構建了二維輪軌滑動溫升有限差分計算模型,分別討論了對流換熱系數(shù)、法向載荷、蠕滑率以及行車速度對鋼軌表面溫升的影響,可為輪軌滑動溫升的準確預測提供參考.

1 計算模型

在列車運行過程中,經過一段時間的磨損和變形,輪軌接觸斑會變得比較狹長,可近似為矩形形狀,從而可在二維框架下討論輪軌滑動溫升問題[1].忽略車輪和鋼軌的內部熱源,考慮熱導率和比熱容的溫度相關性,基于傅里葉導熱定律,二維情況下的輪軌摩擦導熱方程可以寫成

其中,λ(T),c(T) 和ρ 分別表示熱導率、比熱容和密度,溫度場T是位置坐標(x,z) 和時間t的函數(shù),即T=T(x,z;t).令a表示接觸斑縱向寬度的一半,(X,Z)為固定在鋼軌表面的絕對坐標系,(x,z)是原點處于接觸斑前沿的隨動坐標系,如圖1 所示.

圖1 二維輪軌接觸示意圖Fig.1 Schematic of a two-dimensional wheel-rail contact system

當行車速度v保持恒定時,坐標(X,Z)和(x,z)滿足關系式(X,Z)=(x-vt,z),即

將式(2)代入方程(1)中,可得

在穩(wěn)態(tài)情況下,?T/?t=0,方程(3)退化為

對方程(4)中的各物理量進行無量綱處理

在接觸區(qū)域內,輪軌之間因相對滑動而發(fā)生摩擦做功,假設因輪軌滑動摩擦而損失的機械能全部轉化為熱能,并以熱傳導的方式平均傳遞給車輪和鋼軌,忽略輪軌間的相互導熱作用,而接觸區(qū)之外的表面則與空氣發(fā)生對流換熱.因此,邊界條件可表示為

式中,q為熱通量;μ(T)為摩擦系數(shù);h為對流換熱系數(shù);T0為環(huán)境溫度;vs為輪軌間的相對滑動速度,且vs=sv,s為縱向蠕滑率;p(x)表示接觸壓力.

當輪軌材料滿足雙線性各向同性強化本構關系時,接觸壓力p(x)可以寫成[27]

式中,表示Hertz 彈性接觸下使得接觸斑寬度達到2a時所對應的接觸壓力峰值;py為臨界接觸壓力,且py≈1.6σy,σy為屈服強度[28].該模型假設當接觸壓力峰值超過py時,接觸斑內開始出現(xiàn)局部塑性區(qū)域[27].ap為接觸斑內塑性區(qū)域縱向寬度的一半

位于區(qū)域0 ≤|x-a| ＜ap中的表面材料已發(fā)生塑性變形,而位于區(qū)域ap≤|x-a| ≤a中的表面材料仍處于彈性接觸狀態(tài).H*和E*分別滿足

式中,ν1和ν2分別表示車輪和鋼軌的泊松比,E1和E2分別為車輪和鋼軌的彈性模量,H1和H2分別為車輪和鋼軌的塑性硬化模量.

對式(9)在接觸區(qū)域進行積分,可得

式中,R為車輪滾動半徑;Pz為橫向單位長度上的法向載荷,Pz=P/b0,P表示總法向載荷,b0為矩形接觸斑的橫向等效寬度.當H*=E*時,方程(12)退化為

這與Hertz 彈性接觸理論在二維接觸狀態(tài)下的計算結果完全一致[1,29].

由于速度v的影響,車輪高速行進時的總法向載荷(動載)會高于車輪靜止時的總法向載荷(靜載),二者之間滿足[30]

式中,Ps為靜載,等于軸重之半.

2 計算模型的求解

觀察式(7)和式(8)可以發(fā)現(xiàn),同時考慮熱導率、比熱容和摩擦系數(shù)的溫度相關性,使得方程(7)成為具有復雜邊界條件(8) 的含多個變系數(shù)的非線性拋物型偏微分方程,直接進行求解十分困難.因此,基于基爾霍夫變換法,對熱導率的溫度相關性函數(shù)進行積分

從而方程(7)可以轉化成一種更簡潔的形式

其中?-1(T)表示?(T)的反函數(shù).

類似地,邊界條件(8)也可以改寫為

真實材料的κ(?) 和μ(?) 需要通過擬合實驗結果獲得,可能并不總是滿足線性函數(shù)等簡單形式,這意味著攝動法以及拉普拉斯變換等常見方法可能不再適用于方程(16)的求解.因此,具有更強適用性的隱式差分法被用于該問題的求解.

基于以差商近似代替微商的有限差分思想,可將方程(16)和(17)轉化成對應的差分方程形式

圖2 計算區(qū)域的離散Fig.2 Schematic of calculation domain discretization

將方程(19) 代入方程(18) 中,便可以推導出非線性熱傳導方程(16)的二層隱式差分格式

在上述推導過程中,并沒有提前給定溫度相關性函數(shù)λ(T),c(T)和μ(T)的形式,這表明所給出的隱式差分格式(20)可以適用于具有任意溫度相關性函數(shù)形式的熱導率、比熱容和摩擦系數(shù).聯(lián)立求解方程(15)和(20),便可獲得任意節(jié)點位置處的溫度,i=1,2,···,m+1,j=1,2,···,n+1.

需要說明的是,隱式差分法沒有計算區(qū)域維度的限制,當三維彈塑性接觸壓力分布函數(shù)已知時,上述計算方法可較容易地推廣至三維輪軌滑動接觸溫升的求解.

3 結果與討論

3.1 材料屬性的確定

考慮CHN60N 型鋼軌和輪徑(2R) 為890 mm的LMA型車輪.環(huán)境溫度T0=25°C,材料密度ρ=7790 kg/m3,彈性模量E1=E2=209 GPa,硬化模量H1=H2=20.9 GPa,泊松比ν1=ν2=0.3,屈服強度σy=608 MPa.根據(jù)文獻[31]中的實驗數(shù)據(jù),熱導率λ(T)、比熱容c(T)和摩擦系數(shù)μ(T)的擬合曲線如圖3(a)～圖3(c)所示,對應的擬合函數(shù)為

其中,擬合系數(shù)a11=6.61×10-6,a12=0.01,a13=47.71,a21=0.13,a22=487.50,a31=0.35,a32=-1.40×10-3.

圖3 熱導率、比熱容和摩擦系數(shù)的擬合曲線Fig.3 The fitting curves of thermal conductivity,specific heat capacity and friction coefficient

盡管實際材料的屈服強度和彈性模量也會表現(xiàn)出溫度相關性,但其溫度相關性函數(shù)的弓入會大大增加彈塑性接觸問題的復雜性,導致接觸壓力分布函數(shù)的顯式形式不再存在,進而使輪軌滑動溫升的計算變得更加困難.因此,出于簡化問題的考慮,這里忽略了屈服強度和彈性模量的溫度相關性.

3.2 不同模型間的比較

不同輪軌滑動溫升計算模型下的鋼軌表面溫升曲線如圖4 所示.圖中上標(ep),(e),(de)和(in)分別表示彈塑性接觸、彈性接觸、溫度相關性和溫度無關性,后文亦同.由于鋼軌的摩擦熱影響區(qū)深度很淺,僅為2 mm 左右[32-33],為了兼顧計算精度與工作效率,此處可先試取計算區(qū)域深度zmax=4 mm.可以發(fā)現(xiàn),當忽略材料屬性的溫度相關性時,所提出方法計算的溫升曲線與文獻[8]得到的曲線符合較好,表明本文所采用的隱式差分法和選取的計算區(qū)域深度是可行的.同時,考慮材料屬性溫度相關性時的溫升曲線明顯偏離忽略材料屬性溫度相關性時的曲線,二者峰值溫升之間的差值甚至可達640°C.由此可見,在輪軌溫升計算模型中考慮材料屬性的溫度相關性是非常有必要的.另外,在較高溫升下,摩擦系數(shù)會出現(xiàn)劇烈的衰減現(xiàn)象(見圖3(c)),對應的熱通量隨之降低,進而對鋼軌表面溫度的繼續(xù)升高產生抑制作用,最終造成考慮摩擦系數(shù)溫度相關性時的溫升計算結果遠遠低于忽略該特性時的計算結果.因此,在輪軌滑動溫升模擬中十分有必要考慮摩擦系數(shù)的溫度相關性.

圖4 不同輪軌滑動溫升計算模型下的鋼軌表面溫升曲線(Ps=90 kN,h=100 W·(m2·°C)-1,v=350 km/h, s=8%)Fig.4 The temperature rise curves of rail surface for different calculation models of wheel-rail sliding temperature rise(Ps=90 kN,h=100 W·(m2·°C)-1,v=350 km/h, s=8%)

另外,與Hertz 彈性接觸理論相比,當法向載荷過大、接觸斑內表面材料出現(xiàn)局部塑性變形時,彈塑性接觸理論所給出的峰值接觸壓力更小、接觸斑寬度更大(見圖5),進而造成鋼軌表面峰值溫度所在位置朝接觸斑后沿方向移動.由于所采用的高強度鋼軌CHN60N 具有很高的臨界接觸壓力,相應的塑性變形區(qū)域很小,因而兩種接觸壓力形式下的最大滑動溫升比較接近.但對于屈服強度或臨界接觸壓力較小的接觸體材料,這兩種溫升計算結果間的差異將更加顯著.

圖5 彈性接觸壓力與彈塑性接觸壓力的比較Fig.5 The comparison of elastic contact pressure and elasto-plastic contact pressure

3.3 對流換熱系數(shù)對滑動溫升的影響

在輪軌滑動溫升計算模型中,鋼軌非接觸區(qū)域表面的邊界條件直接取決于對流換熱系數(shù),然而該物理參數(shù)往往受到車速、風速以及環(huán)境溫度等諸多因素的影響,其具體取值尚無定論[34].Wang 等[35]通過比例模型實驗測得了鋼軌表面對流換熱系數(shù)在列車高速行駛下的取值范圍為80～160 W/(m2·°C).因此,圖6 和圖7 給出了不同對流換熱系數(shù)下的鋼軌表面溫升和熱通量分布曲線.

圖6 不同對流換熱系數(shù)下的鋼軌表面溫升曲線(Ps=90 kN,v=350 km/h, s=8%)Fig.6 The temperature rise curves of rail surface with different convection coefficients(Ps=90 kN,v=350 km/h, s=8%)

圖7 不同對流換熱系數(shù)下的鋼軌表面熱通量分布曲線(Ps=90 kN,v=350 km/h, s=8%)Fig.7 The heat flux distribution of rail surface with different convection coefficients(Ps=90 kN,v=350 km/h, s=8%)

不難看出,鋼軌表面峰值溫度一直出現(xiàn)在接觸斑內,并靠近接觸斑后沿,這與文獻[36-37]的結論相一致.同時,不同對流換熱系數(shù)下的溫升曲線幾乎完全重合,表明當列車高速行駛時,所討論的對流換熱系數(shù)變化對輪軌滑動溫升的影響甚微,這與文獻[38]的結論一致.進一步觀察圖7 可知:在車輪高速行進情況下,摩擦功率對熱通量的影響會顯著高于對流換熱系數(shù)變化的影響,進而體現(xiàn)為輪軌滑動溫升對對流換熱系數(shù)變化的不敏感性,因此在后續(xù)的討論中對流換熱系數(shù)將被設為定值100 W/(m2·°C).

3.4 法向載荷對滑動溫升峰值的影響

不同法向載荷下的鋼軌表面溫升峰值如圖8 所示.從圖8 中可以發(fā)現(xiàn),在不同計算模型下,對應的溫升峰值均會隨法向載荷的增大近似呈線性上升,但溫升取值則存在明顯的差異.與其他計算模型相比,忽略材料屬性溫度相關性時的溫升峰值計算結果最大,表明忽略材料屬性的溫度相關性會大大高估輪軌間的滑動溫升,考慮材料屬性的溫度相關性可以減小這種高估程度,進而提高溫升計算結果的準確性.同時,僅考慮摩擦系數(shù)溫度相關性時的溫升峰值計算結果明顯低于僅考慮熱導率和比熱容的溫度相關性時的結果,與同時考慮熱導率、比熱容和摩擦系數(shù)的溫度相關性時的溫升計算結果更為接近,表明就式(22) 所表征的材料而言,摩擦系數(shù)的溫度相關性對最終溫升計算結果的影響顯著強于熱導率和比熱容.然而,即便已經在溫升計算模型中弓入了溫度相關的摩擦系數(shù),繼續(xù)考慮熱導率和比熱容的溫度相關性仍會影響最終的溫升峰值計算結果,二者之間的差值可接近70°C,這再一次凸顯了同時考慮多種材料屬性的溫度相關性對于提高輪軌滑動溫升模擬結果準確性的意義.

圖8 鋼軌表面最大溫升隨法向載荷的變化(v=350 km/h, s=8%)Fig.8 The maximum temperature rises of rail surface versus vertical loads(v=350 km/h, s=8%)

3.5 蠕滑率對滑動溫升峰值的影響

鋼軌表面溫升峰值隨蠕滑率的變化曲線如圖9所示.可以看出,溫升峰值隨蠕滑率的增大而上升,這是因為在相同接觸壓力和行車速度的情況下,蠕滑率的增大會弓起輪軌接觸點對間相對速度的增大,提高摩擦熱通量,進而導致鋼軌表面溫度的上升.當蠕滑率較低(s≤1%)時,輪軌滑動溫升并不顯著,熱導率、比熱容和摩擦系數(shù)的變化量很小,不同計算模型的溫升峰值曲線接近重合; 但隨著蠕滑率的進一步增大,鋼軌表面溫度達到較高水平,材料屬性的溫度相關性開始凸顯,因而該4 組曲線互相偏離,且偏離差異隨蠕滑率的增加而增大.忽略材料屬性溫度相關性時的鋼軌滑動溫升峰值與蠕滑率呈線性關系,這與文獻[8] 中的結論相一致; 考慮材料屬性溫度相關性時的溫升峰值-蠕滑率曲線則體現(xiàn)出非線性特征,且位于考慮材料屬性溫度無關性時的曲線的下方.

圖9 鋼軌表面最大溫升隨蠕滑率的變化(Ps=90 kN,v=350 km/h)Fig.9 The maximum temperature rises of rail surface versus creepages(Ps=90 kN,v=350 km/h)

3.6 行車速度對滑動溫升峰值的影響

圖10 鋼軌表面最大溫升隨行車速度的變化(Ps=90 kN, s=8%)Fig.10 The maximum temperature rises of rail surface versus train speeds(Ps=90 kN, s=8%)

圖10 給出了鋼軌表面溫升峰值隨行車速度的變化曲線.與蠕滑率對鋼軌溫升峰值的影響相類似,當法向載荷和蠕滑率恒定時,行車速度的增大會弓起輪軌接觸點對間相對速度的增大,提高摩擦熱通量,進而導致鋼軌表面溫升峰值的增加.因此,熱導率、比熱容和摩擦系數(shù)在高速下的變化量比低速情況更為顯著,從而導致不同溫升模型計算結果之間的差異隨行車速度的增大而加大.在輪軌滑動溫升計算模型中弓入材料屬性的溫度相關性函數(shù),可以避免做出過分高估的溫升評價,這一特征在圖10 中表現(xiàn)為考慮材料屬性溫度相關性時的溫升峰值計算結果明顯低于忽略溫度相關性時的計算結果.同時,僅考慮摩擦系數(shù)溫度相關性時的溫升峰值曲線遠遠低于僅考慮熱導率和比熱容溫度相關性時的結果,并接近于同時考慮這三者的溫度相關性時的曲線,表明圖3(c)所呈現(xiàn)出的摩擦系數(shù)熱衰減特性對溫升計算結果的影響顯著高于熱導率和比熱容.因此,在計算模型中增加摩擦系數(shù)的溫度相關性函數(shù)可大幅提高溫升預測的準確性.

盡管所構建的二維輪軌滑動接觸溫升計算模型同時考慮了接觸壓力的塑性修正和熱導率、比熱容以及摩擦系數(shù)的溫度相關性,但忽略了屈服強度和彈性模量的溫度相關性等因素對溫升結果的影響,與真實的輪軌傳熱狀態(tài)尚有一定的偏差;同時,為了便于獲得接觸壓力的解析解,材料模型采用了簡化的雙線性各向同性強化模型來考慮車輪和鋼軌的彈塑性響應,而實際的車輪和鋼軌的彈塑性響應是非線性的,這將對計算結果產生一定的影響.上述不足之處在目前提出的計算模型中尚無法被考慮,需要在進一步的工作中予以完善.

4 結論

在二維彈塑性接觸模型框架下,同時考慮熱導率、比熱容和摩擦系數(shù)的溫度相關性,基于基爾霍夫變換方法,將復雜的非線性Fourier 導熱方程轉化成一種以熱導率關于溫度的積分函數(shù)為待求量的簡單方程形式,進而構建了相應的輪軌滑動接觸溫升計算模型.基于隱式差分法,推導出了一種不限材料溫度相關性函數(shù)形式的統(tǒng)一差分求解格式,并分別討論了對流換熱系數(shù)、法向載荷、蠕滑率以及行車速度對鋼軌表面滑動溫升的影響,主要結論如下:

(1)由于鋼軌溫升峰值總是位于接觸斑內的后沿附近,并且當接觸斑內出現(xiàn)塑性變形區(qū)域時,彈塑性接觸理論所給出的接觸斑寬度較彈性接觸理論更大,相應的溫升峰值所在位置距接觸斑前沿更遠.

(2)當行車速度處于較高水平時,對流換熱系數(shù)對輪軌滑動溫升的影響甚微,在溫升計算模型中將對流換熱系數(shù)設為定值是可行的.蠕滑率和行車速度的增加,均會弓起輪軌接觸點對間相對速度的增大,摩擦熱通量隨之上升,進而導致鋼軌表面溫升峰值的增加.另外,鋼軌的溫升峰值也會隨法向載荷的增加而近似線性上升.

(3)當實際的輪軌滑動溫升處于較低水平時,材料屬性的溫度相關性對溫升計算結果的影響并不明顯;當實際溫升較顯著時,在計算模型中同時考慮多種材料屬性的溫度相關性可以避免過分高估的溫升評價,且摩擦系數(shù)的溫度相關性對溫升計算結果的影響會明顯高于熱導率和比熱容.