考虑材料温度相关性的二维轮轨弹塑性滑动接触温升分析1)

伏培林 丁立 赵吉中 张旭 阚前华,2) 王平

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

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

引言

近年来,随着我国轨道交通行业的迅速发展,列车的行驶速度和运量均不断提升,这对列车的牵弓和制动性能提出了越来越高的要求.列车的牵弓和制动性能通过轮轨之间的接触关系来实现,一旦列车的牵弓力或制动力超过轮轨接触关系所能提供的粘着力时,轮轨间便会发生相对滑动[1].在滑动过程中,轮轨接触斑内的温度会因摩擦生热而迅速上升,并随着车轮的进一步滑动脱离接触而快速衰减,导致热应力出现交变变化.这种交变的温升诱发热应力,叠加在轴重弓起的机械应力上,将大幅度降低钢轨和车轮的服役寿命[2-3].同时,当滑动温升足够高时,轮轨表层材料发生马氏体相变,并伴随着后续的快速冷却过程而形成脆硬的马氏体组织.在轮轨接触载荷的作用下,这些脆硬组织容易发生脱落,造成轮轨表面的扁疤和材料剥离,进而影响列车行驶的安全性[4-5].另外,过高的轮轨表面温度也会导致材料硬度等力学性能的显著下降,从而对轮轨表面的耐磨损性能带来不利影响[6].由于对轮轨滑动接触温升进行实验测量往往比较困难,因此通过数值模拟和解析分析等非测试方式开展相应的温升研究具有重要的工程意义.

随着ABAQUS、ANASYS 等商业有限元软件的普及,轮轨滑动接触温升的有限元模拟研究受到了学者们的广泛关注.然而,由于轮轨接触斑非常狭小,为了准确地模拟接触斑内的温度场,在有限元模拟过程中必须对这一区域进行大规模网格加密处理,从而所花费的时间成本较为高昂.因此,针对轮轨滑动接触温升问题的解析或半解析研究也得到了大量开展.例如,Tanvir[7]忽略了传热参数的温度相关性,将Hertz 接触压力的椭圆型分布函数近似成四阶多项式级数形式,通过拉普拉斯变换给出了二维轮轨滑动摩擦温升分布的解析解.基于相同的材料属性温度无关性假设,Knothe 和Liebelt[8]通过拉普拉斯变换和格林函数法,给出了任意接触压力下的轮轨滑动温升问题的解析解,讨论了法向接触压力分布形式和表面粗糙度对钢轨表面温度的影响.基于Hertz 弹性接触理论,Ertz 和Knothe[9]讨论了轮轨间换热系数等因素对轮轨摩擦温升的影响.Fischer 等[10]同时考虑了摩擦功和局部塑性变形对接触界面热源的贡献,给出了轮轨滑动接触闪温的解析表达式.此外,考虑到轮轨间存在磨屑等杂物,Wu 等[11]基于三体微接触模型和接触温升理论,讨论了杂物尺寸、杂物密度、表面粗糙度、相对速度以及法向载荷等参数对接触界面温度场的影响.

由于材料属性的温度相关性,摩擦副材料的热力学性质在摩擦温升过程中也会发生变化,考虑材料属性的温度相关性将进一步提升温升计算结果的准确性.然而,材料属性温度相关性的弓入会导致热传导方程的非线性化,大大增加了温升求解的困难性.因此,考虑材料属性的温度相关性,必须首先求解非线性热传导方程.Ling 和Rice[12]考虑了热导率和比热容的温度相关性,通过傅里叶余弦变换,给出了表面具有移动热源半无限大体的准稳态温度场的迭代数值解法.这种求解方法允许任意形式的温度相关性函数,但相应的计算成本高,在实际应用中受到一定的限制.梁钰等[13-14]基于特征正交分解方法,对非线性瞬态热传导问题的有限元离散格式进行降阶处理,从而提出了一种求解考虑传热参数温度相关性的二维非线性瞬态热传导问题的高效算法.通过假设组份材料和孔隙填充物的热导率均随温度呈三次函数变化,戴婷等[15]采用微分求积法求解了一类含孔隙变厚度功能梯度圆盘的稳态温度场.当热扩散系数随温度线性变化时,Fabre 和Hristov[16]基于热平衡积分法给出了一类一维热传导方程的解析解.林府标等[17]通过李群分析法,考虑扩散率关于温度的幂函数形式,得到了该类一维广义热传导方程的显示解析解和行波解.Ghasemi 等[18]则假定热导率是温度的线性函数,将其斜率作为扰动小量,也给出了一类一维非线性热传导问题的摄动解法.此外,微分变换法[19]、有限差分法[20-21]和近场动力学法[22]也被应用于非线性热传导温度场的求解.

由上述分析可见,已有研究存在不足:(1) 接触压力通常被假定为满足Hertz 弹性接触理论中的半椭圆型分布形式,而在实际轮轨接触载荷作用下,接触斑内将发生塑性变形,其接触压力不再满足半椭圆型分布形式,必须进行塑性修正[23-25].(2) 材料参数往往被简化为温度无关的常量,或仅考虑某一个材料参数的温度相关性.实际上,包含摩擦系数在内的多种材料属性均会随温度的改变而发生变化[26],在温升计算模型中同时考虑多种材料属性的温度相关性,有望提升温升计算结果的准确性.

因此,本文通过拟合实验数据,获得了热导率、比热容以及摩擦系数的温度相关性函数,并基于弹塑性接触理论,构建了二维轮轨滑动温升有限差分计算模型,分别讨论了对流换热系数、法向载荷、蠕滑率以及行车速度对钢轨表面温升的影响,可为轮轨滑动温升的准确预测提供参考.

1 计算模型

在列车运行过程中,经过一段时间的磨损和变形,轮轨接触斑会变得比较狭长,可近似为矩形形状,从而可在二维框架下讨论轮轨滑动温升问题[1].忽略车轮和钢轨的内部热源,考虑热导率和比热容的温度相关性,基于傅里叶导热定律,二维情况下的轮轨摩擦导热方程可以写成

其中,λ(T),c(T) 和ρ 分别表示热导率、比热容和密度,温度场T是位置坐标(x,z) 和时间t的函数,即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)中,可得

在稳态情况下,∂T/∂t=0,方程(3)退化为

对方程(4)中的各物理量进行无量纲处理

在接触区域内,轮轨之间因相对滑动而发生摩擦做功,假设因轮轨滑动摩擦而损失的机械能全部转化为热能,并以热传导的方式平均传递给车轮和钢轨,忽略轮轨间的相互导热作用,而接触区之外的表面则与空气发生对流换热.因此,边界条件可表示为

式中,q为热通量;μ(T)为摩擦系数;h为对流换热系数;T0为环境温度;vs为轮轨间的相对滑动速度,且vs=sv,s为纵向蠕滑率;p(x)表示接触压力.

当轮轨材料满足双线性各向同性强化本构关系时,接触压力p(x)可以写成[27]

式中,表示Hertz 弹性接触下使得接触斑宽度达到2a时所对应的接触压力峰值;py为临界接触压力,且py≈1.6σy,σy为屈服强度[28].该模型假设当接触压力峰值超过py时,接触斑内开始出现局部塑性区域[27].ap为接触斑内塑性区域纵向宽度的一半

位于区域0 ≤|x-a| ＜ap中的表面材料已发生塑性变形,而位于区域ap≤|x-a| ≤a中的表面材料仍处于弹性接触状态.H*和E*分别满足

式中,ν1和ν2分别表示车轮和钢轨的泊松比,E1和E2分别为车轮和钢轨的弹性模量,H1和H2分别为车轮和钢轨的塑性硬化模量.

对式(9)在接触区域进行积分,可得

式中,R为车轮滚动半径;Pz为横向单位长度上的法向载荷,Pz=P/b0,P表示总法向载荷,b0为矩形接触斑的横向等效宽度.当H*=E*时,方程(12)退化为

这与Hertz 弹性接触理论在二维接触状态下的计算结果完全一致[1,29].

由于速度v的影响,车轮高速行进时的总法向载荷(动载)会高于车轮静止时的总法向载荷(静载),二者之间满足[30]

式中,Ps为静载,等于轴重之半.

2 计算模型的求解

观察式(7)和式(8)可以发现,同时考虑热导率、比热容和摩擦系数的温度相关性,使得方程(7)成为具有复杂边界条件(8) 的含多个变系数的非线性抛物型偏微分方程,直接进行求解十分困难.因此,基于基尔霍夫变换法,对热导率的温度相关性函数进行积分

从而方程(7)可以转化成一种更简洁的形式

其中ϑ-1(T)表示ϑ(T)的反函数.

类似地,边界条件(8)也可以改写为

真实材料的κ(ϑ) 和μ(ϑ) 需要通过拟合实验结果获得,可能并不总是满足线性函数等简单形式,这意味着摄动法以及拉普拉斯变换等常见方法可能不再适用于方程(16)的求解.因此,具有更强适用性的隐式差分法被用于该问题的求解.

基于以差商近似代替微商的有限差分思想,可将方程(16)和(17)转化成对应的差分方程形式

图2 计算区域的离散Fig.2 Schematic of calculation domain discretization

将方程(19) 代入方程(18) 中,便可以推导出非线性热传导方程(16)的二层隐式差分格式

在上述推导过程中,并没有提前给定温度相关性函数λ(T),c(T)和μ(T)的形式,这表明所给出的隐式差分格式(20)可以适用于具有任意温度相关性函数形式的热导率、比热容和摩擦系数.联立求解方程(15)和(20),便可获得任意节点位置处的温度,i=1,2,···,m+1,j=1,2,···,n+1.

需要说明的是,隐式差分法没有计算区域维度的限制,当三维弹塑性接触压力分布函数已知时,上述计算方法可较容易地推广至三维轮轨滑动接触温升的求解.

3 结果与讨论

3.1 材料属性的确定

考虑CHN60N 型钢轨和轮径(2R) 为890 mm的LMA型车轮.环境温度T0=25°C,材料密度ρ=7790 kg/m3,弹性模量E1=E2=209 GPa,硬化模量H1=H2=20.9 GPa,泊松比ν1=ν2=0.3,屈服强度σy=608 MPa.根据文献[31]中的实验数据,热导率λ(T)、比热容c(T)和摩擦系数μ(T)的拟合曲线如图3(a)～图3(c)所示,对应的拟合函数为

其中,拟合系数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 热导率、比热容和摩擦系数的拟合曲线Fig.3 The fitting curves of thermal conductivity,specific heat capacity and friction coefficient

尽管实际材料的屈服强度和弹性模量也会表现出温度相关性,但其温度相关性函数的弓入会大大增加弹塑性接触问题的复杂性,导致接触压力分布函数的显式形式不再存在,进而使轮轨滑动温升的计算变得更加困难.因此,出于简化问题的考虑,这里忽略了屈服强度和弹性模量的温度相关性.

3.2 不同模型间的比较

不同轮轨滑动温升计算模型下的钢轨表面温升曲线如图4 所示.图中上标(ep),(e),(de)和(in)分别表示弹塑性接触、弹性接触、温度相关性和温度无关性,后文亦同.由于钢轨的摩擦热影响区深度很浅,仅为2 mm 左右[32-33],为了兼顾计算精度与工作效率,此处可先试取计算区域深度zmax=4 mm.可以发现,当忽略材料属性的温度相关性时,所提出方法计算的温升曲线与文献[8]得到的曲线符合较好,表明本文所采用的隐式差分法和选取的计算区域深度是可行的.同时,考虑材料属性温度相关性时的温升曲线明显偏离忽略材料属性温度相关性时的曲线,二者峰值温升之间的差值甚至可达640°C.由此可见,在轮轨温升计算模型中考虑材料属性的温度相关性是非常有必要的.另外,在较高温升下,摩擦系数会出现剧烈的衰减现象(见图3(c)),对应的热通量随之降低,进而对钢轨表面温度的继续升高产生抑制作用,最终造成考虑摩擦系数温度相关性时的温升计算结果远远低于忽略该特性时的计算结果.因此,在轮轨滑动温升模拟中十分有必要考虑摩擦系数的温度相关性.

图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 弹性接触理论相比,当法向载荷过大、接触斑内表面材料出现局部塑性变形时,弹塑性接触理论所给出的峰值接触压力更小、接触斑宽度更大(见图5),进而造成钢轨表面峰值温度所在位置朝接触斑后沿方向移动.由于所采用的高强度钢轨CHN60N 具有很高的临界接触压力,相应的塑性变形区域很小,因而两种接触压力形式下的最大滑动温升比较接近.但对于屈服强度或临界接触压力较小的接触体材料,这两种温升计算结果间的差异将更加显著.

图5 弹性接触压力与弹塑性接触压力的比较Fig.5 The comparison of elastic contact pressure and elasto-plastic contact pressure

3.3 对流换热系数对滑动温升的影响

在轮轨滑动温升计算模型中,钢轨非接触区域表面的边界条件直接取决于对流换热系数,然而该物理参数往往受到车速、风速以及环境温度等诸多因素的影响,其具体取值尚无定论[34].Wang 等[35]通过比例模型实验测得了钢轨表面对流换热系数在列车高速行驶下的取值范围为80～160 W/(m2·°C).因此,图6 和图7 给出了不同对流换热系数下的钢轨表面温升和热通量分布曲线.

图6 不同对流换热系数下的钢轨表面温升曲线(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 不同对流换热系数下的钢轨表面热通量分布曲线(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%)

不难看出,钢轨表面峰值温度一直出现在接触斑内,并靠近接触斑后沿,这与文献[36-37]的结论相一致.同时,不同对流换热系数下的温升曲线几乎完全重合,表明当列车高速行驶时,所讨论的对流换热系数变化对轮轨滑动温升的影响甚微,这与文献[38]的结论一致.进一步观察图7 可知:在车轮高速行进情况下,摩擦功率对热通量的影响会显著高于对流换热系数变化的影响,进而体现为轮轨滑动温升对对流换热系数变化的不敏感性,因此在后续的讨论中对流换热系数将被设为定值100 W/(m2·°C).

3.4 法向载荷对滑动温升峰值的影响

不同法向载荷下的钢轨表面温升峰值如图8 所示.从图8 中可以发现,在不同计算模型下,对应的温升峰值均会随法向载荷的增大近似呈线性上升,但温升取值则存在明显的差异.与其他计算模型相比,忽略材料属性温度相关性时的温升峰值计算结果最大,表明忽略材料属性的温度相关性会大大高估轮轨间的滑动温升,考虑材料属性的温度相关性可以减小这种高估程度,进而提高温升计算结果的准确性.同时,仅考虑摩擦系数温度相关性时的温升峰值计算结果明显低于仅考虑热导率和比热容的温度相关性时的结果,与同时考虑热导率、比热容和摩擦系数的温度相关性时的温升计算结果更为接近,表明就式(22) 所表征的材料而言,摩擦系数的温度相关性对最终温升计算结果的影响显著强于热导率和比热容.然而,即便已经在温升计算模型中弓入了温度相关的摩擦系数,继续考虑热导率和比热容的温度相关性仍会影响最终的温升峰值计算结果,二者之间的差值可接近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%)时,轮轨滑动温升并不显著,热导率、比热容和摩擦系数的变化量很小,不同计算模型的温升峰值曲线接近重合; 但随着蠕滑率的进一步增大,钢轨表面温度达到较高水平,材料属性的温度相关性开始凸显,因而该4 组曲线互相偏离,且偏离差异随蠕滑率的增加而增大.忽略材料属性温度相关性时的钢轨滑动温升峰值与蠕滑率呈线性关系,这与文献[8] 中的结论相一致; 考虑材料属性温度相关性时的温升峰值-蠕滑率曲线则体现出非线性特征,且位于考虑材料属性温度无关性时的曲线的下方.

图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 给出了钢轨表面温升峰值随行车速度的变化曲线.与蠕滑率对钢轨温升峰值的影响相类似,当法向载荷和蠕滑率恒定时,行车速度的增大会弓起轮轨接触点对间相对速度的增大,提高摩擦热通量,进而导致钢轨表面温升峰值的增加.因此,热导率、比热容和摩擦系数在高速下的变化量比低速情况更为显著,从而导致不同温升模型计算结果之间的差异随行车速度的增大而加大.在轮轨滑动温升计算模型中弓入材料属性的温度相关性函数,可以避免做出过分高估的温升评价,这一特征在图10 中表现为考虑材料属性温度相关性时的温升峰值计算结果明显低于忽略温度相关性时的计算结果.同时,仅考虑摩擦系数温度相关性时的温升峰值曲线远远低于仅考虑热导率和比热容温度相关性时的结果,并接近于同时考虑这三者的温度相关性时的曲线,表明图3(c)所呈现出的摩擦系数热衰减特性对温升计算结果的影响显著高于热导率和比热容.因此,在计算模型中增加摩擦系数的温度相关性函数可大幅提高温升预测的准确性.

尽管所构建的二维轮轨滑动接触温升计算模型同时考虑了接触压力的塑性修正和热导率、比热容以及摩擦系数的温度相关性,但忽略了屈服强度和弹性模量的温度相关性等因素对温升结果的影响,与真实的轮轨传热状态尚有一定的偏差;同时,为了便于获得接触压力的解析解,材料模型采用了简化的双线性各向同性强化模型来考虑车轮和钢轨的弹塑性响应,而实际的车轮和钢轨的弹塑性响应是非线性的,这将对计算结果产生一定的影响.上述不足之处在目前提出的计算模型中尚无法被考虑,需要在进一步的工作中予以完善.

4 结论

在二维弹塑性接触模型框架下,同时考虑热导率、比热容和摩擦系数的温度相关性,基于基尔霍夫变换方法,将复杂的非线性Fourier 导热方程转化成一种以热导率关于温度的积分函数为待求量的简单方程形式,进而构建了相应的轮轨滑动接触温升计算模型.基于隐式差分法,推导出了一种不限材料温度相关性函数形式的统一差分求解格式,并分别讨论了对流换热系数、法向载荷、蠕滑率以及行车速度对钢轨表面滑动温升的影响,主要结论如下:

(1)由于钢轨温升峰值总是位于接触斑内的后沿附近,并且当接触斑内出现塑性变形区域时,弹塑性接触理论所给出的接触斑宽度较弹性接触理论更大,相应的温升峰值所在位置距接触斑前沿更远.

(2)当行车速度处于较高水平时,对流换热系数对轮轨滑动温升的影响甚微,在温升计算模型中将对流换热系数设为定值是可行的.蠕滑率和行车速度的增加,均会弓起轮轨接触点对间相对速度的增大,摩擦热通量随之上升,进而导致钢轨表面温升峰值的增加.另外,钢轨的温升峰值也会随法向载荷的增加而近似线性上升.

(3)当实际的轮轨滑动温升处于较低水平时,材料属性的温度相关性对温升计算结果的影响并不明显;当实际温升较显著时,在计算模型中同时考虑多种材料属性的温度相关性可以避免过分高估的温升评价,且摩擦系数的温度相关性对温升计算结果的影响会明显高于热导率和比热容.