2019, Vol. 47
2. 上海市轨道交通结构耐久与系统安全重点实验室，上海 201804;
3. 比亚迪汽车工业有限公司，广东 深圳 518118;
4. 朔黄铁路发展有限责任公司，河北 沧州 062356
2. Shanghai Key Laboratory of Rail Infrastructure Durability and System Safety, Shanghai 201804, China;
3. BYD Automotive Industry Co. Ltd., Shenzhen 518118, China;
4. Shuohuang Railway Development Co. Ltd., Cangzhou 062356, China
重载列车速度和轴重的不断提高，加剧了重载铁路钢轨的磨耗.磨耗使得轮轨接触状态发生改变，轮轨接触由单点接触变为两点接触或共形接触的概率大大提升，轮轨接触界面的作用变得更加复杂.研究表明，钢轨磨耗与轮轨法向接触应力直接相关[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理论中最小余能方程为
| $ \begin{array}{l} \min C = \frac{1}{2}{p_{Iz}}{A_{IzJz}}{p_{Jz}} + \left( {{h_I} - \delta } \right){p_{Iz}}\\ \quad {\rm{s}}{\rm{.t}}{\rm{.}}\;\;{p_{Jz}} \ge 0, {S_0}\sum\limits_{J = 1}^M {{p_{Jz}}} = {F_z} \end{array} $ | (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单元的影响系数，处理过程遵循以下步骤：
|
图 2 影响系数计算 Fig.2 Calculation of influencing coefficient |
(1) 拟合钢轨型面，选取可能接触区域内I单元，确定其法向方向nI和切向方向sI；
(2) 对作用在J单元的法向单位力F按照nI和sI分解为FnI和FsI；
(3) 设I单元和J单元的切线斜率分别为kI和kJ，以确定I单元和J单元的法向夹角α.
完成(1)~(3)步骤后，可以近似认为I单元在某一个平面内受到一个法向力FnI和一个切向力FsI的作用，由于是在局部坐标系中的计算，故局部坐标系的影响系数可以表示为
| $ {A_{InJn}} = {A_{IzJz}}\cos \alpha - {A_{IzJy}}\sin \alpha $ | (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) = {b_1}{x^2} + {b_2}xy + {b_3}{y^2} + {b_4}x + {b_5}y + {b_6} $ | (3) |
式中：(x，y)为单元的全局坐标；b1、b2、…、b6均为根据车轮和钢轨型面确定的参数.
这种方法仍然是基于弹性半空间假设，其认为钢轨接触表面为平面或者近似平面，轮轨在轨距角接触时显然是不成立的.本文采用在局部坐标系中的一种近似处理的方法，如图 4所示，将车轮与钢轨视作二维平面内的两条曲线，不考虑纵向的影响，计算可能接触区域内P单元的法向间隙，处理过程遵循以下步骤：
|
图 4 法向间隙的计算 Fig.4 Calculation of normal distance |
(1) 拟合车轮与钢轨型面，对钢轨表面进行网格划分，选取可能接触区域内的P单元，作P单元的法线，交车轮型面曲线于Q0点；
(2) 选取角度θ，将直线PQ0以P为旋转中心逆时针旋转θ，交车轮型面曲线于Q1点，计算PQ1；在PQ0另一侧的车轮型面曲线上寻求点Q′1，使得PQ1=PQ′1；
(3) 作∠Q1PQ′1的平分线，交车轮型面曲线于Q点，计算PQ，若|PQ－PQ0|≤Δ，Δ表示精度，则认为PQ0的长度即为要求的法向间隙；若不满足，则将Q点视作Q0点，重复(1)~(3)步骤，直到满足条件.
为验证这种方法的准确性，采用1.2节中的计算参数，补充的参数为选取LMA系列原始车轮型面，轮对横移量为+12 mm，θ取0.1°，Δ取10-5 mm，利用有限元方法和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 |
当车轮的横移量为0时，轮轨接触为轨顶面接触，属于平面接触范畴，计算结果如图 7所示.此时轮轨接触斑为近似椭圆形，接触面积为196.83 mm2，最大法向应力为1 757 MPa.
|
图 7 横移量为0时的轮轨法向接触应力 Fig.7 Normal contact stress at 0 mm lateral displacement |
当轮对横移量为+13 mm时，轮轨接触为轨距角接触，此时为两点接触，计算结果如图 8所示.此时接触斑形状不再是椭圆形，两个接触斑的面积分别为69.66 mm2和165.24 mm2，对应的最大法向接触应力分别为1 582 MPa和1 740 MPa.
|
图 8 横移量为+13 mm时的轮轨法向接触应力 Fig.8 Normal contact stress at +13 mm lateral displacement |
由于列车运行时是蛇形运动，因此轮对存在一定的横移量，尤其在曲线运行时，更是如此.利用迹线法确定在此钢轨磨耗状态下不同横移量时的轮轨接触点，选取-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 |
|
图 13 靠近轨道中心线的接触斑最大法向应力和接触面积 Fig.13 Maximum contact stress and contact area of the region closer to the track center line |
基于Kalker法向接触的最小余能方程，对影响系数和法向间隙进行修正，引入法向角计算非平面接触时的影响系数，在局部坐标系下获得接触点的法向间隙，所采用的方法与有限元方法得到的结果接近，验证了方法的准确性.对修正后的最小余能方程进行编程，可以求解轮轨接触的应力分布，并选取某重载铁路通过总重达100 Mt的CHN75型面磨耗钢轨，车轮选取LMA系列原始型面，利用所编程序，计算了30 t轴重下重载磨耗钢轨的轮轨法向接触特性，得到以下结论：
(1) 当横移量为-12~+12 mm时，轮轨接触为单点接触，轮轨接触状态接近，最大法向应力和接触面积小幅变化，且最大法向应力与接触面积的变化趋势相反.
(2) 当横移量为+12.9~+13.2 mm时，轮轨接触状态为两点接触，在两个接触斑上最大法向应力和接触面积变化幅度较大，且最大法向应力与接触面积的变化趋势相同.
(3) 当横移量达到+14 mm时，轮轨接触为车轮轮缘和钢轨轨距角接触，此时接触斑变得狭长，接触面积较小，但最大接触应力达到3 287 MPa.
本文对钢轨磨耗状态下的法向接触特性进行了研究，但轮轨法向接触应力与钢轨磨耗之间的关系还有待更深入探讨，这是笔者未来研究的重点.同时，如何让计算结果更适用于车辆动力学仿真也是需要进一步研究的方向.
| [1] |
侯传伦, 翟婉明, 邓锐. 曲线磨耗状态下轮轨弹塑性接触有限元分析[J]. 中国铁道科学, 2009, 30(5): 28 HOU Chuanlun, ZHAI Wanming, DENG Rui. Finite element analysis of the elastic-plastic contact of the worn wheels and rails on curve[J]. China Railway Science, 2009, 30(5): 28 DOI:10.3321/j.issn:1001-4632.2009.05.005 |
| [2] |
王平, 刘哲. 轮轨磨耗对岔区轮轨接触行为的影响[J]. 铁道工程学报, 2016, 33(1): 46 WANG Ping, LIU Zhe. Effects of wheel and rail wear on wheel/rail contact behavior in turnout zone[J]. Journal of Railway Engineering Society, 2016, 33(1): 46 DOI:10.3969/j.issn.1006-2106.2016.01.010 |
| [3] |
KALKER J J. Three-dimensional elastic bodies in rolling contact[M]. Dordrecht: Kluwer Academic Publishers, 1990
|
| [4] |
KALKER J J. A fast algorithm for the simplified theory of rolling contact[J]. Vehicle System Dynamics, 1982, 11(1): 1 DOI:10.1080/00423118208968684 |
| [5] |
KALKER J J. A simplified theory for non-hertzian contact[J]. Vehicle System Dynamics, 1983, 12(1/2/3): 43 |
| [6] |
LI Z. Wheel-rail rolling contact and its application to wear simulation[D]. Delft: Delft University of Technology, 2002.
|
| [7] |
BURGELMAN N, LI Z, DOLLEVOET R. A new rolling contact method applied to conformal contact and the train-turnout interaction[J]. Wear, 2014, 321: 94 DOI:10.1016/j.wear.2014.10.008 |
| [8] |
金学松.轮轨蠕滑理论及其试验研究[D].成都: 西南交通大学, 1999. JIN Xuesong. The theory of wheel rail creepage and its experimental study[D]. Chengdu: Southwest Jiaotong University, 1999. |
| [9] |
REN Z, IWNICKI S D, XIE G. A new method for determining wheel-rail multi-point contact[J]. Vehicle System Dynamics, 2011, 49(10): 1533 |
| [10] |
杨新文, 顾少杰, 练松良. 30 t轴重重载列车轮轨法向接触应力分析[J]. 铁道学报, 2015(6): 19 YANG Xinwen, GU Shaojie, LIAN Songliang. An analysis of wheel-rail normal contact stress under 30 t-axle-load heavy haul train[J]. Journal of the China Railway Society, 2015(6): 19 DOI:10.3969/j.issn.1001-8360.2015.06.003 |
| [11] |
VOLLEBREGT E, SEGAL G. Solving conformal wheel-rail rolling contact problems[J]. Vehicle System Dynamics, 2014, 52(S1): 455 |
| [12] |
MA Y, MARKINE V, MASHAL A A, et al. Modelling verification and influence of operational patterns on tribological behaviour of wheel-rail interaction[J]. Tribology International, 2017, 114: 264 DOI:10.1016/j.triboint.2017.04.038 |
| [13] |
孙宇, 翟婉明. 钢轨磨耗演变预测模型研究[J]. 铁道学报, 2017, 39(8): 1 SUN Yu, ZHAI Wanming. A prediction model for rail wear evolution[J]. Journal of the China Railway Society, 2017, 39(8): 1 DOI:10.3969/j.issn.1001-8360.2017.08.001 |
| [14] |
杨新文, 姚一鸣, 周顺华. 基于轮轨非Hertz接触的影响系数的有限元计算方法[J]. 同济大学学报(自然科学版), 2017, 45(10): 1476 YANG Xinwen, YAO Yiming, ZHOU Shunhua. Calculation of influencing number of wheel-rail non-Hertz contact using finite element method[J]. Journal of Tongji University(Natural Science), 2017, 45(10): 1476 DOI:10.11908/j.issn.0253-374x.2017.10.009 |
| [15] |
BLANCO-LORENZO J, SANTAMARIA J, VADILLO E G, et al. On the influence of conformity on wheel-rail rolling contact mechanics[J]. Tribology International, 2016, 103: 647 DOI:10.1016/j.triboint.2016.07.017 |