复杂运营条件下重载货车车轮磨耗发展的数值预测

2019-04-01王树国

同济大学学报（自然科学版） 2019年1期

王璞，王树国

(中国铁道科学研究院集团有限公司铁道建筑研究所，北京 100081)

铁路运输中，列车车轮磨耗现象一直是最为普遍且难以解决的一个问题.磨耗引起车轮型面的改变，会导致轮轨接触状态和系统动力性能发生变化，对列车运行稳定性、安全性以及车辆轨道系统各部件的使用寿命都会产生较大影响.对重载铁路而言，轮轨磨耗问题尤为突出，养修工作巨大，磨耗超过一定限度如未及时镟修或更换，则可能引起重大事故.

通过数值仿真对重载列车车轮磨耗发展进行预测，可指导车辆部门有针对性地安排和优化养护维修工作，实现从发现问题后进行治理性养护到预判问题开展预防性养护的转变，有效减小养修工作量和成本，减缓车轮磨耗发展.另外，通过磨耗发展定量预测，可作为车辆系统和轨道结构设计的评价标准，为轮轨型面优化、车辆轨道结构改进提供科学指导.因此，轮轨磨耗数值预测具有重要意义，国内外学者在这方面进行了一系列的研究工作[1-12].然而，既有研究对于轨道条件的考虑通常过于简化，实际运营中列车可能会以多种速度在不同轨道条件地段运行，既有方法对这种复杂的运行条件考虑很少，亦或只是近似反映轨道的统计特征，与实际存在较大差异.另外，预测方法多采用特定运行距离作为车轮型面更新步长，即列车运行固定距离后车轮型面进行更新.这种方法一定程度上可提高计算效率，但可能会在计算中引起磨耗的偏差不断累积、放大甚至直接导致错误发生，很难确保数值计算的可靠性和稳定性.

本文建立了重载列车车轮磨耗发展的数值预测模型，并基于MATLAB编制了计算程序.模型由车辆-轨道耦合动力学计算、轮轨滚动接触分析、材料摩擦损失模型、复杂运营条件的模拟策略、型面更新自适应步长方法等模块组成.基于Archard材料磨损理论，在车辆轨道耦合动力学仿真和轮轨局部接触分析基础上进行轮轨接触斑磨耗深度分布计算；通过多工况仿真并引入量纲为一权重因子来实现对重载铁路实际复杂运营条件的合理模拟；采用了一种自适应步长算法进行车轮型面更新，可有效减少误差累积、改善数值模型稳定性和可靠性.最后，利用本文模型对大秦重载铁路实际运营条件下列车车轮磨耗发展过程进行数值预测分析.

1 重载车辆-轨道耦合动力学模型

车轮磨耗发展预测是在车辆-轨道耦合动力学和轮轨滚动接触计算基础上进行的，首先建立重载车辆-轨道耦合动力学模型.

基于多体动力学理论，建立三大件式转向架重载货车车辆模型.对车体及转向架结构部件(摇枕、侧架、轴箱、轮对等)均采用6自由度刚体模拟，对车体与心盘、旁承，楔块减振器与摇枕、侧架，以及侧架与轴箱间的接触摩擦作用均采用点-面接触单元组合进行模拟，并充分考虑部件实际间隙.如对心盘与车体间的接触摩擦作用模拟如图1所示，车体垂向荷载导致的法向接触力以及由于摇枕车体相对转动产生的摩擦力通过均匀分布于心盘表面的8个点-面接触单元进行模拟，当车体与心盘部分出现脱离产生间隙时，相应位置接触摩擦作用为0.中央悬挂枕簧采用双线性刚度阻尼力学单元模拟，通过两个分别对应内簧和外簧的刚度矩阵和内外簧高度差参数来描述其特殊的非线性特性.通过精细化建模确保模型尽可能与实际情况相符.三大件式转向架模型如图2所示.

轮轨接触计算基于Hertz接触理论[13]及Kalker的FASTSIM算法[14]进行，主要分为接触点位置探测、整体接触力学量计算(接触力、蠕滑率等)以及接触斑局部接触力学量计算(接触应力、蠕滑应力、滑动速度等).

轮轨接触计算中考虑两点接触的影响，当车轮沿着钢轨横向移动过程中，轮轨接触状况变化如图3所示.最初轮轨在走行面上发生一点接触，随着车轮横移，在某一临界横移量时发生走行面和轮缘两点接触，当车轮继续横移时，出现爬轨现象，只在轮缘上发生一点接触，走行面范围内轮轨脱离.值得注意的是，车轮连续横移过程中，两点接触的情形只出现在某一临界位置，而且出现前后接触点位置发生突变，因此可以根据此进行两点接触位置的探测.

图1 心盘-车辆接触摩擦作用Fig.1 Center plate-vehicle contact friction interaction

图2 三大件式转向架力学模型Fig.2 Mechanical model of three-piece bogie

图3 轮轨接触点变化过程Fig.3 Changing process of wheel-rail contact point

对于车轮某一特定侧滚角，逐渐改变车轮横移量，分别计算接触点位置，当接触点位置发生明显突变时，表明两点接触发生在上一次车轮横移改变量范围内，将此范围再细分为多个离散段，进而寻找两点接触发生的更小的区间范围，然后再进行细分……重复几次上述细分过程，可以获得足够精确的该侧滚角条件下两点接触发生的位置.同样的方法可得到所有侧滚角条件下的两点接触位置及对应的车轮横移量.

轮轨接触作用力与下部轨道结构密不可分，模型充分考虑轨道的弹性和阻尼特性，刚度阻尼参数基于道床垂向刚度和横向阻力测试以及扣件刚度阻力测试结果[15-17]，最后得到整体轨道刚度阻尼参数取值.轨道力学模型如图4所示，N1、N2分别为走行面接触和轮缘接触位置的法向接触力，T1、T2分别为两接触位置的蠕滑力，Gy、Gz分别为轨道系统对钢轨的垂横向约束反力，Cy、Cz、Dy和Dz分别为轨道系统的垂横向刚度和阻尼.

a 单点接触b 两点接触

图4轨道力学模型

Fig.4Mechanicalmodeloftrack

轮轨接触计算中，为了求得法向力和蠕滑力分布，将接触斑进行离散化处理，如图5所示.接触斑坐标系原点位于接触斑中心，x、y轴分别对应轮轨接触面的纵、横向.首先将接触斑沿y轴等分成ny个长条，然后将每个长条均等分成nx个矩形离散块，则每个矩形块尺寸为

(1)

式中：(x,y)为离散块中心在接触斑坐标系下的位置坐标；a、b分别为接触斑半长轴、半短轴长度；dy、dx(y)分别为离散块沿y方向、x方向的长度，dx(y)由离散块y坐标决定.这种划分方法可提高接触斑边界网格分辨率和结果精确性.接触斑经以上离散化后，通过迭代计算求得各网格上的蠕滑力及蠕滑速度，即得接触斑蠕滑力和蠕滑速度分布，同时完成接触斑滑动区、黏着区划分.通过将各网格内的蠕滑力累加得接触斑总蠕滑力.

图5 轮轨接触斑离散化Fig.5 Discretization of wheel-rail contact patch

2 车轮磨耗计算模型

基于车辆-轨道耦合动力学模型算得轮轨接触斑几何特征、接触斑滑动区和黏着区分布以及接触法向力和蠕滑速度分布计算接触斑内的磨耗深度分布.接触斑内网格划分与轮轨局部接触计算时保持一致.

对接触斑每个离散单元，基于Archard材料磨损理论[18]计算其磨耗体积损失，如图6所示，ΔVw为离散单元上的磨耗体积损失量，ΔDw为磨耗深度，v为离散单元处轮轨表面相对滑动速度.

ΔVw(x,y)=kw(x,y)

(2)

式中：(x,y)为离散单元中心在接触斑坐标系下的位置坐标；Δs(x,y)为离散单元处对应的轮轨相对滑动距离；kw(x,y)为磨耗系数，根据离散单元处法向应力和相对滑动速度由图7得到[2].p(x,y)为离散单元中心位置的法向应力，根据赫兹接触理论，计算如下：

(3)

式中，N为接触斑上总的法向接触力.

对于轮轨相对滑动距离Δs(x,y)，在轮轨接触计算过程中可得到接触斑内黏着区、滑动区分布，如果所考虑的离散单元位于黏着区，则滑动距离Δs(x,y)=0，无磨耗发生；如果离散单元位于滑动区，Δs(x,y)计算如下：

(4)

式中：V0为车轮滚动速度.离散单元处轮轨相对滑动速度v(x,y)计算如下：

v(x,y)=vr(x,y)-ve(x,y)=

(5)

式中：vr(x,y)表示刚体相对滑动速度分量，ve(x,y)为弹性变形相对滑动速度分量；ξx、ξy、φ为接触斑纵、横向、自旋蠕滑率，u(x,y)=(ux(x,y),uy(x,y))为离散单元处弹性变形位移量.通常弹性变形滑动分量ve(x,y)比刚体滑动分量vr(x,y)小得多，故计算中将其忽略.

图6 轮轨接触斑磨耗深度分布计算模型Fig.6 Calculation model of wear distribution in wheel-rail contact patch

图7 磨耗系数取值Fig.7 Value of wear coefficient

通过以上计算，可得到接触斑每个离散单元上的磨耗体积损失ΔVw(x,y)，进而由下式可得到每个离散单元对应的磨耗深度ΔDw(x,y)，即接触斑内磨耗深度分布如下：

(6)

将具有相同y坐标的离散单元处的磨耗深度相加，可得到接触斑内每一纵向长条上对应的磨耗深度，即为磨耗深度沿轮轨型面的分布.

对动力学仿真中的每一积分步，通过以上计算过程均能得到对应的磨耗深度沿轮轨型面的分布结果(接触斑坐标系下)，然后结合当前积分步的轮轨接触点位置信息(接触点在车轮坐标系下的水平位置坐标，动力学仿真可得)将磨耗分布定位到车轮型面坐标系下，物理意义为车轮当前位置型面在车轮滚动一周时所产生的磨耗(车轮滚动一周，当前位置型面与钢轨只接触一次).磨耗计算中单次动力学仿真均设置为一种车辆以固定速度在单一轨道条件下运行，因此各积分步对应的运行条件是相同的，所引起的对应位置车轮磨耗情况也应是相似的.因此，对所考虑时间范围内各积分步的磨耗结果求平均，结果意义为当前运行条件下车轮每滚动一周所引起的型面磨耗深度分布，为方便应用将其转换为车轮滚动单位距离所引起的磨耗深度分布，定义为车轮磨耗速率cw，如图8所示，cw为车轮型面坐标系水平坐标yw的函数.

图8 车轮磨耗速率Fig.8 Wheel wear rate

3 重载铁路复杂运营条件的模拟

重载铁路的实际运营条件较为复杂，列车可能会以不同的速度在由多种轨道条件地段组成的铁路线路上运行，这都会对车轮磨耗发展情况产生显著影响，并大大增加了数值预测的难度.为对这种复杂运营条件进行模拟，根据列车行驶过程中可能存在的轨道条件和运行速度设计工况，对每种工况分别进行动力学仿真和磨耗计算，得到各车轮的磨耗速率(车轮滚动单位距离所产生的型面磨耗深度分布)，然后引入速度权重因子α和轨道条件权重因子β确定各工况所占比率，最后得到列车在实际运营条件下各车轮的加权平均磨耗速率.

假定所考虑铁路线路存在m种轨道条件(T1,T2,T3,…,Tm)；在其中任一Tj轨道条件地段，列车可能有nj种运行速度(V1j,V2j,V3j,…,Vnjj)；各车轮的加权平均磨耗速率计算如下：

(7)

式中：下标k为轮对编号，l、r分别表示左、右侧车轮；cw,k,l,ij(yw)、cw,k,r,ij(yw)分别表示列车在Tj轨道条件地段以Vij速度运行工况下各车轮的磨耗速率；βw,j为Tj轨道条件对应的权重因子，αw,ij为Tj轨道条件下列车运行速度Vij对应的权重因子.βw,j取值由所考虑铁路线路状况决定，有

(8)

式中：S(Tj)表示Tj轨道条件地段总长度，与线路总长度比值为βw,j值.αw,ij根据行车组织资料及实际调研统计确定，为列车在Tj轨道条件地段以速度Vij运行的频率.

4 车轮磨耗发展的型面更新策略

轮轨动力相互作用是车轮磨耗产生的根源，磨耗引起车轮踏面材料损失会造成车轮型面的改变，反过来又会对轮轨接触几何关系和动力相互作用产生显著影响.因此车轮磨耗发展是一个交互作用的过程，车轮型面逐渐连续变化.而在数值模拟时是无法对连续过程进行仿真的，必须进行离散化处理.因此将车轮磨耗发展考虑为离散的过程，通过迭代计算实现，在每一迭代步内假设车轮型面固定不变，型面变化引起的轮轨动力行为变化忽略不计，磨耗发展简化为线性变化.每一迭代步最后，由磨耗速率及步长信息计算车轮累积磨耗，更新车轮型面，新型面作为下一迭代步计算的初始型面输入.

计算过程中迭代步步长选择是一个重要问题.步长意义为一个迭代步内车辆运行距离.步长过小会造成巨大计算代价，而步长过大将降低数值模型的精确和稳定性.固定步长可使数值计算更简单高效，但如果在某一迭代步计算中出现意外偏差算得的磨耗速率过大，由于步长固定，相应将会造成当前步累积磨耗过大，导致错误的计算结果或直接导致后续计算无法进行.

因此本文采用一种自适应步长的车轮型面更新算法，设定每一迭代步车轮累积磨耗峰值达到固定限值时，迭代步终止，执行型面更新.各迭代步的步长根据磨耗速率而不断调整变化，这种自适应步长策略可有效避免错误的出现，同时能够减小误差累积，改善数值模型的可靠性和稳定性.下面具体介绍这种型面更新算法，每一迭代步的具体计算过程如下：

首先，通过动力学仿真和磨耗计算得到实际运营条件下车辆各个车轮的加权平均磨耗速率cw,k,l(yw)、cw,k,r(yw)，然后计算所有车轮的磨耗速率峰值为

(9)

根据迭代步内磨耗线性发展假设，在每一迭代步中车辆各个车轮的磨耗速率是固定不变的.设置各迭代步累积磨耗深度限值为ξw，则迭代步内列车运行距离(迭代步步长)计算如下：

P=ξw/cw,max

(10)

在确定步长信息后，迭代步内每个车轮的累积磨耗可计算如下：

(11)

式中：Cw,k,l(yw)、Cw,k,r(yw)表示迭代步内各个车轮的累积磨耗深度分布，为车轮型面坐标系水平坐标yw的函数.

最后，根据累积磨耗Cw,k,l(yw)、Cw,k,r(yw)对车轮型面进行更新，更新后的型面作为初始型面代入下一迭代步进行计算.值得注意的是，Cw,k,l(yw)、Cw,k,r(yw)通常情况下并非平滑的分布曲线，而是存在尖锐形状变化(主要由高非线性轮轨滚动接触计算引起)，进而造成更新后的型面也对应出现这种尖锐变化，与现场实际情况存在差异.因此，每一迭代步有必要在车轮型面更新后进行适当的平滑处理.本文采用B样条函数对更新后的型面进行平滑处理，然后将处理后的型面代入下一迭代步计算.

综上，可建立起完整的重载列车车轮磨耗发展数值预测方法，其涉及车辆-轨道耦合动力学仿真、轮轨局部接触分析、接触斑及车轮型面磨耗分布计算、重载铁路复杂运营条件模拟策略以及自适应步长的车轮型面更新算法等内容.图9清晰地呈现了整个计算过程.

5 大秦铁路实际运营条件下重载列车车轮磨耗发展预测分析

大秦铁路是我国第一条双线电气化开行重载单元列车的运煤专线，全线桥梁、隧道众多，并存在各种不同半径曲线.基于本文模型对大秦铁路(湖东站—柳村站)实际运营条件下列车车轮的磨耗发展过程进行数值预测.

首先对线路分布及行车速度情况进行统计，如表1所示.分别对不同运行条件工况进行动力学仿真和磨耗计算，然后基于第3节方法对实际运营条件进行模拟.大秦铁路上、下行线路有所区别，但下行线运行的空车车轮踏面磨耗远远小于重车，因此只考虑上行线路重车车轮磨耗发展情况，空车返回时产生的磨耗忽略不计.

图9 车轮磨耗发展数值预测的迭代计算流程Fig.9 Iterative computation process for numerical prediction of wheel wear development

值得注意的是，实际线路中曲线地段相对于列车行进方向有左曲线也有右曲线，在对相应工况进行计算时，统一按照左曲线进行计算，然后将同轮对左右轮算得的磨耗速率平均处理，作为该轮对车轮的磨耗速率指标.具体做法如下：

列车在Tj曲线轨道条件地段以Vij速度运行工况条件下，计算得到的第k轮对左、右车轮的磨耗速率为cw,k,l,ij(yw)、cw,k,r,ij(yw)，采用式(12)进行修正如下：

cw,k,l′,ij(yw)=cw,k,r′,ij(yw)=

(12)

然后基于修正后的各车轮磨耗速率cw,k,l′,ij(yw)、cw,k,r′,ij(yw)，通过权重因子αw,ij和βw,j计算列车实际运行时各车轮的加权平均磨耗速率如下：

(13)

这是一种简化处理方法，对线路条件进行调研时统计明确各种半径的左曲线和右曲线分布，分别建立动力学模型进行计算是更为精确的方法，但通常会耗费较大计算代价.

下面对C80型重车车辆在大秦铁路上运行时各个车轮的磨耗发展过程进行预测分析.为方便表述，对各个车轮进行编号，如图10所示.根据大秦铁路实际情况，车轮型面采用LM磨耗型踏面，钢轨为75 kg·m-1钢轨，轨底坡1/40.采用自适应步长的车轮型面更新策略，设置每一迭代步终止的累积磨耗深度限值为0.05 mm.预测计算结果如图11和图12所示，其中图11给出了不同运行里程条件下各车轮的磨耗发展分布情况及磨耗导致的型面变化，图12给出了各车轮最大累积磨耗深度的增长曲线.

表1 大秦铁路线路分布及行车速度(总长620 109 m)Tab.1 The line distribution and vehicle speed of Da-Qin railway (overall length 620 109 m)

图10 车轮编号Fig.10 Wheel number

可以看出，随着车辆运行里程的增加，车辆各个车轮的磨耗均不断增大，但磨耗发展呈逐渐减缓的趋势.各个车轮的磨耗主要分布在名义滚动圆两侧走行区域，起导向作用的1、2号车轮(1位轮对)和5、6号车轮(3位轮对)磨耗分布范围相比3、4、7、8号车轮的磨耗分布范围更宽.各车轮在靠近轮缘侧的磨耗发展均更快，对导向轮对车轮(1、2、5、6号)这一特征更为明显，1、2、5、6号车轮的最大磨耗深度位置更靠近轮缘而3、4、7、8号车轮的最大磨耗深度位置更靠近名义滚动圆.具体如下：

1、2号车轮磨耗主要分布在[-35 mm，40 mm]范围，随着运行里程的增加，最大磨耗深度均主要发生在[20 mm，25 mm]范围内，当运行里程达到35 000 km时，最大磨耗深度达2.197 mm.