重载铁路钢轨磨耗状态下的轮轨法向接触特性

许玉德，严道斌，孙小辉，唐永康

(1.同济大学 道路与交通工程教育部重点实验室，上海 201804；2.上海市轨道交通结构耐久与系统安全重点实验室，上海 201804；3.比亚迪汽车工业有限公司，广东 深圳 518118；4.朔黄铁路发展有限责任公司，河北 沧州 062356)

重载列车速度和轴重的不断提高，加剧了重载铁路钢轨的磨耗.磨耗使得轮轨接触状态发生改变，轮轨接触由单点接触变为两点接触或共形接触的概率大大提升，轮轨接触界面的作用变得更加复杂.研究表明，钢轨磨耗与轮轨法向接触应力直接相关[1-2].

国内外学者对轮轨接触理论展开了深入研究.Kalker[3-5]提出三维弹性体滚动接触理论，基于此开发的CONTACT程序，可以求解考虑蠕滑的轮轨接触问题；为了适用车辆动力学仿真，Kalker提出简化理论并开发了FASTSIM程序.Li[6-7]对Kalker在CONTACT中采用的弹性半空间影响系数作了修正，提出四分之一空间的影响系数有限元计算方法，同时在计算轮轨接触时不再将蠕滑视作常数，基于此开发的WEAR程序，能够有效处理轨距角轮轨共形接触问题.金学松[8]在Kalker理论的基础上，分析了轮轨滚动接触蠕滑率与蠕滑力之间的关系，拓展了其简化理论模型，使其更适合求解轮轨黏着问题.任尊松[9]在迹线法的基础上提出了迹线极值法，建立车辆-轨道系统动力学模型，能够解决轮轨多点接触问题并获得准确的轮轨多点接触几何参数.杨新文[10]提出基于赫兹理论的P_M(partition-model)模型，采用切片投影法寻求接触区域，并利用分区方法计算非椭圆接触斑上的法向接触应力.

当钢轨产生磨耗时，车轮与钢轨的接触可能会发生在轨距角区域，尤其是在小半径曲线上.轨距角区域的轮轨接触是曲面接触，弹性半空间假设不再成立[11].CONTACT程序由于采用弹性半空间假设，因此不能较好地处理这类问题；而有限元方法不受弹性半空间假设的限制，计算结果准确，是目前处理轨距角接触问题常用的选择[12]，其缺点是计算效率较低.本文基于Kalker三维弹性体滚动接触理论，引入法向角对弹性半空间假设影响系数进行修正，在局部坐标系下计算轮轨接触的法向间隙，对修正后的Kalker法向接触方程编程求解，可以准确、快速地计算轮轨法向接触应力，提出的方法适用于轮轨平面和非平面接触.

1 轮轨法向接触应力计算

1.1 轮轨滚动接触理论

Kalker的三维弹性体滚动接触理论，是一种基于边界元的数值计算方法.如图1所示，该方法将可能的接触区域划分为矩形单元网格，利用弹性半空间Boussinesq-Cerruti公式计算得到任一单元与可能接触区域内其他所有单元的位移差，借助Newton-Raphson方法对接触方程迭代求解得到所有接触单元的应力分布[13].

图1 可能接触区域的网格划分Fig.1 Meshing of possible contact region

对于法向接触问题，Kalker理论中最小余能方程为

(1)

式中：C为余能；I、J为矩形单元的编号；pIz为作用在I单元的法向应力；AIzJz为影响系数，表示作用在J单元的法向单位力在I单元处引起的法向位移；hI为I单元处轮轨的法向间隙；δ为法向侵入量；S0为每个矩形单元的面积，S0=Δx1·Δx2；M为接触区域内单元数量；Fz为法向力.

1.2 影响系数的修正

影响系数AIzJz对法向接触方程的求解精度起到重要作用，在法向接触问题研究中影响系数的计算多采用Boussinesq-Cerruti弹性半空间解析解.当钢轨产生磨耗时，车轮与钢轨可能在轨距角区域接触，弹性半空间假设不再成立.有限元计算影响系数的方法不受弹性半空间假设限制，因此计算影响系数的结果是准确的.文献[14]利用有限元方法计算影响系数，但计算效率较低.文献[6]认为弹性半空间假设在处理法向接触问题时仍然是一种合理的选择.文献[15]提出了一种考虑法向角的非平面影响系数计算方法.

如图2所示，将钢轨视作二维平面内的曲线，不考虑纵向方向的影响，计算J单元受法向单位力作用下可能接触区域内I单元的影响系数，处理过程遵循以下步骤：

(1)拟合钢轨型面，选取可能接触区域内I单元，确定其法向方向nI和切向方向sI；

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

(3)设I单元和J单元的切线斜率分别为kI和kJ，以确定I单元和J单元的法向夹角α.

图2 影响系数计算Fig.2 Calculation of influencing coefficient

完成(1)～(3)步骤后，可以近似认为I单元在某一个平面内受到一个法向力FnI和一个切向力FsI的作用，由于是在局部坐标系中的计算，故局部坐标系的影响系数可以表示为

AInJn=AIzJzcosα-AIzJysinα

(2)

式中：AInJn为局部坐标系的影响系数；AIzJz和AIzJy为全局坐标系的影响系数，由Boussinesq-Cerruti公式计算得到；α为法向角.

以接触区域的曲率半径15 mm为例，进行影响系数计算.选取CHN75钢轨型面，在轨距角处划定一个10 mm×10 mm的矩形可能接触区域，每个单元是0.5 mm×0.5 mm的正方形单元，在矩形区域的中心点处作用一个单位法向力，计算可能接触区域内各单元的影响系数.为验证该近似方法的准确性，利用有限元方法和CONTACT计算在相同条件下的影响系数，结果如图3所示.

图3 影响系数计算方法对比Fig.3 Comparison of calculation methods of influencing coefficient

从影响系数的最大值来看，有限元方法计算所得影响系数最大值为3.06×10-6mm，本文方法为3.10×10-6mm，误差为+1.31%；CONTACT基于Boussinesq-Cerruti弹性半空间解析解的结果为2.92×10-6mm，误差为-4.6%.从变化趋势来看，3种方法计算的结果趋势比较一致，但本文方法与有限元方法得到的结果更接近，因此本文采用的方法是准确的.

1.3 法向间隙的修正

法向间隙是计算法向接触问题的重要输入参数.常用的轮轨法向间隙计算方法是将轮轨接触看作是全局坐标下的接触，然后计算两点之间沿y轴方向的差值，作为法向间隙.CONTACT程序中认为任一单元的法向间隙可用抛物线形函数表征

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

(3)

式中：(x，y)为单元的全局坐标；b1、b2、…、b6均为根据车轮和钢轨型面确定的参数.

这种方法仍然是基于弹性半空间假设，其认为钢轨接触表面为平面或者近似平面，轮轨在轨距角接触时显然是不成立的.本文采用在局部坐标系中的一种近似处理的方法，如图4所示，将车轮与钢轨视作二维平面内的两条曲线，不考虑纵向的影响，计算可能接触区域内P单元的法向间隙，处理过程遵循以下步骤：

(1)拟合车轮与钢轨型面，对钢轨表面进行网格划分，选取可能接触区域内的P单元，作P单元的法线，交车轮型面曲线于Q0点；

图4 法向间隙的计算Fig.4 Calculation of normal distance

为验证这种方法的准确性，采用1.2节中的计算参数，补充的参数为选取LMA系列原始车轮型面，轮对横移量为+12 mm，θ取0.1°，Δ取10-5mm，利用有限元方法和CONTACT计算在相同条件下的法向间隙，结果如图5所示.

图5 法向间隙计算结果对比Fig.5 Comparison of normal distance results

本文方法与有限元方法的计算结果比较接近，尤其是在接触点中心附近区域.从利用CONTACT抛物线多项式差值计算的结果来看，这种方法在处理非平面问题时存在较大误差，尤其是在离接触点中心较远的地方.在横向方向-5 mm处有限元方法法向间隙的结果为0.032 mm，本文方法计算结果为0.031 mm，CONTACT计算结果则为0.051 mm，相比于有限元结果误差达59.4%；在+5 mm处有限元方法法向间隙的结果为0.050 mm，本文方法计算结果为0.049 mm，CONTACT计算结果则为0.071 mm，相比于有限元结果误差达42.0%.本文方法与有限元方法结果较为一致，因此准确性较好.

基于上述影响系数和法向间隙的修正，对Kalker三维弹性体滚动接触理论中法向接触问题的最小余能方程进行编程求解，可以准确、快速地计算轮轨接触的法向应力分布情况.

2 钢轨磨耗状态下的法向接触特性

以30 t轴重列车为例进行计算，车轮选取LMA系列原始型面，钢轨选取CHN75轨通过总重达100 Mt的磨耗型面，如图6所示，该磨耗型面测量于我国某条重载铁路.需要说明的是，本文后续图中，均以钢轨“工”字形截面的水平向为横向方向，铅垂向为垂向方向，沿钢轨长度方向为纵向方向.

图6 重载铁路CHN75磨耗钢轨Fig.6 Heavy haul CHN75 worn rail

2.1 轨顶面接触

当车轮的横移量为0时，轮轨接触为轨顶面接触，属于平面接触范畴，计算结果如图7所示.此时轮轨接触斑为近似椭圆形，接触面积为196.83 mm2，最大法向应力为1 757 MPa.

2.2 轨距角接触

当轮对横移量为+13 mm时，轮轨接触为轨距角接触，此时为两点接触，计算结果如图8所示.此时接触斑形状不再是椭圆形，两个接触斑的面积分别为69.66 mm2和165.24 mm2，对应的最大法向接触应力分别为1 582 MPa和1 740 MPa.

图7 横移量为0时的轮轨法向接触应力Fig.7 Normal contact stress at 0 mm lateral displacement

图8 横移量为+13 mm时的轮轨法向接触应力Fig.8 Normal contact stress at +13 mm lateral displacement

2.3 不同横移量下的法向接触特性

由于列车运行时是蛇形运动，因此轮对存在一定的横移量，尤其在曲线运行时，更是如此.利用迹线法确定在此钢轨磨耗状态下不同横移量时的轮轨接触点，选取-12～+14 mm的横移量，每隔2 mm取一种工况，计算得到不同横移量下的轮轨接触应力和接触面积，见图9.可以看出，在横移量为-12～+10 mm时，接触斑形态变化不大；在+12～+14 mm时，轮轨接触由车轮踏面与钢轨轨顶面接触变为车轮轮缘和钢轨轨距角的接触，最大接触应力达到3 287 MPa.

图9 不同横移量下的法向接触应力Fig.9 Normal contact stress at each lateral displacement

最大法向应力和接触面积的计算结果如图10所示.最大法向应力的变化趋势与接触面积的变化趋势相反，随着接触面积的减小，最大法向应力相应增大.

图10 不同横移量下最大法向应力与接触面积Fig.10 Maximum contact stress and contact area at each lateral displacement

在横移量为+12.9～+13.2 mm区间时出现了两点接触的现象，如图11所示.

图11 不同横移量下两点接触的法向接触应力Fig.11 Normal contact stress at two-point contact

对两点接触的两个接触斑进行分析，最大法向应力和接触面积的结果如图12和图13所示.在两点接触的情况下，最大法向应力与接触面积变化的趋势是一致的.

图12 靠近钢轨中心线的接触斑最大法向应力和接触面积Fig.12 Maximum contact stress and contact area of the region closer to the rail center line

3 结论

基于Kalker法向接触的最小余能方程，对影响系数和法向间隙进行修正，引入法向角计算非平面接触时的影响系数，在局部坐标系下获得接触点的法向间隙，所采用的方法与有限元方法得到的结果接近，验证了方法的准确性.对修正后的最小余能方程进行编程，可以求解轮轨接触的应力分布，并选取某重载铁路通过总重达100 Mt的CHN75型面磨耗钢轨，车轮选取LMA系列原始型面，利用所编程序，计算了30 t轴重下重载磨耗钢轨的轮轨法向接触特性，得到以下结论：

图13 靠近轨道中心线的接触斑最大法向应力和接触面积Fig.13 Maximum contact stress and contact area of the region closer to the track center line

(1)当横移量为-12～+12 mm时，轮轨接触为单点接触，轮轨接触状态接近，最大法向应力和接触面积小幅变化，且最大法向应力与接触面积的变化趋势相反.

(2)当横移量为+12.9～+13.2 mm时，轮轨接触状态为两点接触，在两个接触斑上最大法向应力和接触面积变化幅度较大，且最大法向应力与接触面积的变化趋势相同.

(3)当横移量达到+14 mm时，轮轨接触为车轮轮缘和钢轨轨距角接触，此时接触斑变得狭长，接触面积较小，但最大接触应力达到3 287 MPa.

本文对钢轨磨耗状态下的法向接触特性进行了研究，但轮轨法向接触应力与钢轨磨耗之间的关系还有待更深入探讨，这是笔者未来研究的重点.同时，如何让计算结果更适用于车辆动力学仿真也是需要进一步研究的方向.