道岔作为线路的重要设备，由转辙器、辙叉与连接部件3部分组成. 车辆通过转辙器时，车轮在基本轨与尖轨间转换，会引起剧烈的轮轨相互作用，造成轮轨磨耗. 这会降低列车运行舒适性，缩短钢轨使用寿命甚至危害机车车辆运行安全^[1]. 轮轨滚动接触模型是研究轮轨磨耗的一个关键，轮轨滚动接触模型大致可分为Hertz接触模型与非Hertz接触模型.

Hertz接触理论以及基于Hertz理论的切向接触模型可简单快速求解轮轨接触问题，因此被广泛应用于轮轨动力学的法向、切向力求解. 文献[2-4]均建立了车辆-轨道耦合动力学模型，基于Hertz理论求解轮轨法向力，基于切向算法（FASTSIM）求解轮轨切向力. Hertz接触模型假设接触斑内曲率恒定，接触斑内法向间隙沿接触斑纵横向均对称分布. 但随着轮轨磨耗，出现轮缘与轨距角接触，不能满足Hertz假设，此时采用Hertz模型得到的接触斑与应力分布均有较大误差，因此国内外研究学者提出了非Hertz模型.

几种基于虚拟渗透法的快速非Hertz滚动接触算法被提出，它们假设接触信息在车轮前进方向仍满足赫兹假设，但在接触斑横向不再满足赫兹假设. 文献[5]采用K-P、LINDER、STRIPES、ANALYN非赫兹快速算法研究区间线路轮轨接触斑内磨耗分布， 但以上几种非赫兹快速算法忽略摇头角的影响，因此Liu等^[6-7]在KP基础上分别发展了EKP算法和MKP算法，考虑摇头角后其法向间隙沿接触斑纵横向均为非对称，且大摇头角对接触斑内接触应力分布影响较大. 为更精确地模拟轮轨间接触问题，Kalker^[8]将接触区进行离散，基于弹性半空间假设，利用最小余能求解轮轨间接触问题. 由于Kalker三维非赫兹滚动接触理论精度较高，可考虑任意接触几何，因此广泛被应用于轮轨损伤分析中. 文献[9]利用Kalker三维非赫兹接触理论求解轮轨法向应力分布，分析重载铁路曲线超高对钢轨磨耗的影响. 文献[10]利用Kalker三维非赫兹接触理论计算高速铁路磨耗车轮与60N钢轨静态接触性能，并与CHN60钢轨进行对比，发现60N钢轨在相同条件下，能够减少轮轨滚动接触疲劳损伤的发生. 通过有限元建立的轨道模型虽然可精确模拟轮轨接触几何^[11]，但计算效率太低，不适合大规模的轮轨接触计算.

综上分析，目前采用的轮轨滚动接触解析算法在计算岔区轮轨接触时，未能同时考虑道岔尖轨变截面与摇头角的影响. 因此，本文提出的一种适用于道岔区的三维非对称接触几何算法可考虑两者的影响，将得到的三维轮轨法向间隙代入Kalker三维非赫兹滚动接触理论可提高岔区轮轨接触应力求解的精确性. 此外，本文将三维非对称接触几何算法结合车辆-道岔耦合动力学模型和英国谢菲尔大学提出的磨耗模型（USFD磨耗模型）求解车轮通过曲尖轨时曲尖轨的磨耗分布，并与不考虑摇头角与变截面的几何接触算法得到的结果对比，分析尖轨变截面与摇头角对曲尖轨磨耗的影响，本文的研究结果可为道岔轮轨接触与损伤分析提供参考.

1.   接触几何计算方法

为了分析尖轨变截面与摇头角对轮轨滚动接触行为的影响，本文采用两种算法计算轮轨接触几何：第一种是未考虑尖轨变截面与摇头角的简化算法；第二种是考虑这两种因素的精确算法.

1.1   简化算法

简化算法利用迹线法求解车轮与道岔二维截面的接触点，迹线法得到接触点的详细求解过程可参考文献[12].

根据迹线法得到接触点后，建立接触点局部坐标系${O_{\text{c}}} {\text{-}} {X_{\text{c}}}{Y_{\text{c}}}{Z_{\text{c}}}$，如图1所示，原点位于轮轨接触点，O_cX_c指向车轮滚动方向，O_cZ_c指向接触点法向，O_cY_c由右手法则确定. 将轮轨型面转换到接触点坐标系，求解沿O_cY_c轴的轮轨几何间隙$f({Y_{\text{c}}})$，简化算法假设法向间隙沿接触斑纵轴O_cX_c满足Hertz假设，轮轨三维法向间隙h可根据式（1）得到.

图  1  接触点局部坐标

Figure  1.  Wheel-rail contact coordinate system

下载: 全尺寸图片 幻灯片

式中：R_c为接触点处车轮滚动圆半径.

1.2   精确算法

简化算法采用迹线法求解轮轨接触点，然而迹线法基于等截面钢轨推导得到，且简化算法忽略了摇头角对法向间隙的影响. 对于道岔尖轨，钢轨截面的宽度与高度沿轨道纵向均发生变化，为更贴近真实的接触情况，本文提出精确算法，可考虑摇头角与尖轨变截面几何廓形对轮轨接触行为的影响.

首先，建立轨道坐标系O-XYZ，原点位于轨道中心，OX沿轨道纵向并指向轮对运动方向，OZ指向垂直于轨道平面，OY为横向；其次，建立轮对坐标系O_w-X_wY_wZ_w，O_w固结于轮对质心，随轮对一起运动，O_wX_w指向轮对前进方向，O_wY_w沿轮轴方向，O_wZ_w垂直于O_wX_wY_w平面. 将车轮视为回转体，用X_w= 0的平面切割车轮，其与车轮型面的交线即为车轮主轮廓线，其余轮廓线（X_wi）的横坐标Y_wi与主轮廓线保持一致，则垂坐标为

当轮对摇头角为$\psi$，侧滚角为$\phi$，横移量为$\Delta y$时，车轮轮廓线在轨道坐标系中的表达式为

为获得三维的尖轨型面，本文采用三次样条函数对尖轨关键断面进行插值得到任意顶宽的尖轨型面^[13].

采用三维的轮轨型面，计算量较大. 为了提高计算效率，先采用迹线法得到粗略的接触点位置. 由于尖轨宽度与高度沿轨道方向变化率较小，因此，假设真实接触点在迹线法得到的接触点附近. 本文在得到的粗略接触点附近一定范围内利用车轮三维型面在变截面尖轨上插值求解最小距离从而找到接触点，大大缩短了计算时间. 如图2所示，${d_{{\text{st}}}}$和${d_{{\text{sw}}}}$分别为车轮轮廓线到基本轨和尖轨的最短距离，车轮上点$P$(${X'_{{\text{w}}i}},{Y'_{{\text{w}}i}},{Z'_{{\text{w}}i}}$)对应钢轨上点$Q$(${X'_{{\text{w}}i}},{Y'_{{\text{w}}i}},{Z'_{{\text{r}}i}}$)，其中${Z'_{{\text{r}}i}}$可通过插值得到，如式（4）.

图  2  精确算法求解轮轨接触点示意

Figure  2.  Schematic diagram of the exact method for determining the wheel-rail contact point

下载: 全尺寸图片 幻灯片

式中：${f_{{\text{interpolate}}}}$为插值函数；(${X_{\text{r}}},{Y_{\text{r}}},{Z_{\text{r}}}$)为尖轨变截面上的离散点.

为了提高计算效率，可先确定X =${X'_{{\text{w}}i}}$平面对应的钢轨廓形，利用${Y'_{{\text{w}}i}}$在此钢轨廓形上插值，即可得到钢轨对应的${Z'_{{\text{r}}i}}$.

已知轮轨坐标可求解粗略接触点附近区域内左、右侧每条车轮轮廓线与钢轨的最小距离$\Delta {Z_{{\rm{L}}n\min }}$和$\Delta {Z_{{\rm{R}}n\min }}$，根据式（5）得到左轮与右轮对应的最小距离分别为$\Delta {Z}_{{\rm{L}}\mathrm{min}}$、$\Delta {Z}_{{\rm{R}}\mathrm{min}}$，$\Delta {Z}_{{\rm{L}}\mathrm{min}}\leqslant \Delta {Z}_{{\rm{L}}n\mathrm{min}}$，$\Delta {Z}_{{\rm{R}}\mathrm{min}}\leqslant \Delta {Z}_{{\rm{R}}n\mathrm{min}}$.

得到最小距离后判断车轮是否与左、右轨同时发生接触，判断条件为$\left| {\Delta {Z_{{\rm{Lmin}} }} - \Delta {Z_{{\rm{Rmin}} }}} \right| < {\text{1}}{{\text{0}}^{ - 6}}\;{\text{mm}}$，若不满足判断条件则进行侧滚角迭代，迭代侧滚角为

式中：${Y_{{\rm{L}}{\text{cmin}}}} 、{Y_{{\rm{R}}{\text{cmin}}}}$分别是左、右轮轨接触点的横坐标.

道岔转辙器部件内，轮轨易发生两点接触，尖轨与基本轨会同时承担车轮荷载. 参考文献[14]，可将组合廓形分为尖轨区和基本轨区. 得到主接触点后，进一步得到${d_{{\text{st}}}}$和${d_{{\text{sw}}}}$，式（5）得到的最小距离$\Delta {Z_{{\rm{L}}\min }}$可能是${d_{{\text{st}}}}$或者${d_{{\text{sw}}}}$. 根据轮轨间弹性压缩量确定车轮、尖轨和基本轨是否同时发生接触.

车轮与两股钢轨接触状态为

式中：$\delta$为Hertz渗透量.

确定轮轨接触点后，根据接触点坐标与轮轨坐标的相对距离可将轨道坐标系下的轮轨型面坐标转换到局部坐标系，如式（8）.

式中：(X_cw,r，Y_cw,r，Z_cw,r)为轮轨型面坐标转换到接触点局部坐标系后的坐标；$\Delta {X_{{\text{w}},{\text{r}}}}、\Delta {Y_{{\text{w}},{\text{r}}}}、\Delta {Z_{{\text{w}},{\text{r}}}}$分别为轮轨坐标相对于轮轨接触点的距离.

则道岔区三维法向间隙为

式中：Z_cw、Z_cr分别为车轮和钢轨坐标在接触点坐标系下的Z坐标.

2.   轮轨滚动接触行为分析

利用上文的两种轮轨接触几何算法可得到轮轨法向间隙，然后，通过Kalker非赫兹滚动接触理论得到轮轨接触解，其中，轮轨三维法向间隙是Kalker接触理论的重要输入参数. 以S1002CN车轮踏面与12号道岔曲尖轨为例，车轮名义滚动圆半径为460 mm，轨距为1 435 mm，轨底坡为1/40. 轮轨接触特性的计算参数：轮轨摩擦系数为0.4，泊松比为0.28，剪切模量为82 GPa.

2.1   接触点分析

为了对比两种算法得到的接触点差异，摇头角$\psi$设为0和10 mrad，横移量∆Y设为 −2.0 ~ 8.0 mm，步长为1 mm，尖轨顶宽设为49 mm，变截面顶宽为48.3 ~ 49.7 mm，接触点结果如图3所示.

图  3  不同横移量下的接触点位置

Figure  3.  Positions of the contact points under different lateral displacements

下载: 全尺寸图片 幻灯片

由图3（a）可知：当摇头角为0或10 mrad时，不同横移量下两种算法得到的接触点横坐标差异较小，可忽略不计. 由图3（b）可知：当摇头角为0时，简化算法得到的接触点纵坐标均为0，而精确算法却随着横移量的增大纵坐标呈增大趋势. 这是由于摇头角为0时，轮对无超前和滞后距离，当采用迹线法时，迹线纵坐标始终为0；精确算法考虑了尖轨变截面，尖轨沿轨道纵向高度与宽度变化，因此，最小距离所在的点不一定位于主轮廓线. 当摇头角为10 mrad时，两种算法得到的接触点纵坐标差异较大. 这是由于精确算法考虑尖轨变截面与车轮真实型面，车轮接触的钢轨截面随接触点纵向距离变化而变化；而简化算法假设钢轨为等截面，其接触的钢轨截面始终相同.

2.2   接触斑形状分析

为分析轮对摇头角与尖轨变截面对轮轨接触斑的影响，利用两种算法计算横移量为 −2.0 mm和7.5 mm，$\psi$ = 10 mrad情况下的轮轨接触斑，结果如图4所示.

图  4  接触斑轮廓对比

Figure  4.  Comparison of contact patch contour

下载: 全尺寸图片 幻灯片

由图4 （a）可知：当横移量为 −2.0 mm时，两种算法所得结果差异较小，这是由于车轮踏面与尖轨轨顶接触，尖轨在轨顶处截面变化不大. 由图4 （b）可知：当横移量为7.5 mm时，轮缘与尖轨轨距角接触，两种算法得到的接触斑形状的差异较明显. 精确算法得到的接触斑形状沿接触斑纵横向均为非对称，简化算法得到的接触斑形状关于Y_c轴对称，关于X_c轴非对称；简化算法得到的接触斑面积为76 mm²，精确算法所得结果为83 mm²，比简化算法得到的结果大9.2%. 这是由于尖轨廓形在轨距角位置变化较大，精确算法考虑尖轨变截面后，接触点附近简化算法与精确算法轮轨廓形差异较大，此外，精确算法考虑摇头角影响后，接触斑关于X_c轴有一定偏移.

2.3   接触应力分析

进一步分析摇头角为10 mrad，横移量为7.5 mm时法向接触应力的分布情况，结果如图5所示.

图  5  横移量为7.5 mm时接触斑内法向应力

Figure  5.  Pressure distribution within contact patch at a lateral displacement of 7.5 mm

下载: 全尺寸图片 幻灯片

由图5（a）可知：当X_c= −2 mm时，随着Y_c增大，法向应力差异变大，这是由于偏离接触点的横向距离越大，摇头角和变截面影响变大，两种算法得到的三维法向间隙差异增大. 由图5（b）可知：当Y_c= 8 mm时，简化算法得到的法向应力关于X_c轴对称分布，分布的范围为 −3 ~ 3 mm，其最大法向应力为439 MPa；精确算法结果关于X_c轴非对称分布，分布范围为−6.4 ~ 1.2 mm，其最大法向应力为552 MPa，比简化算法大25.7%. 这是由于两种算法得到的接触斑形状差异较大，因此法向应力分布也差异明显.

3.   轮轨磨耗分析

由第2节分析，轮轨接触几何会显著影响轮轨接触应力与接触斑的计算，而这些是轮轨磨耗分析的重要输入参量. 为了分析尖轨变截面与摇头角对尖轨磨耗的影响，本文首先建立车辆道岔耦合动力学模型仿真分析，得到轮对运动参数以及轮轨间法向力，然后利用两种接触几何算法以及Kalker三维非赫兹滚动接触理论计算相应的轮轨接触解，最后将这些轮轨接触参数代入USFD磨耗模型^[15]分析钢轨磨耗深度.

3.1   车辆道岔耦合动力学模型

以CRH3型动车组与12号道岔曲尖轨为研究对象，采用SIMPACK建立多体动力学模型，包括1个车体、2个转向架、4个轮对和8个轴箱和转臂，其中模型各部件质量、转动惯量以及一系、二系悬挂参数均采用设计参数进行设置. 车辆侧逆向过岔，速度为40 km/h.

3.2   磨耗预测模型

本文采用USFD磨耗预测模型分析接触斑内的磨耗分布. 接触斑内的摩擦功W(X,Y)可以通过磨耗指数表示，如式（10）.

式中：p_X和p_Y分别为接触斑单元的纵向和横向蠕滑力；${\gamma _X}$和${\gamma _Y}$分别为对应的蠕滑率.

磨耗深度可通过式（11）计算得到.

式中：$\rho$为钢轨密度；$\Delta X$为接触斑单元宽度；${K_{\rm{w}}}$为磨耗系数，如式（12）.

参考文献[16]，利用快速叠加法求解某一钢轨截面上的磨耗量. 假定车轮滚动通过道岔某断面过程中接触斑信息不变，因此道岔某断面的磨耗量可根据此接触斑磨耗分布沿线路方向叠加得到. 假定接触斑沿线路纵向和横向划分的单元格个数分别为N₁和N₂，则钢轨截面沿横向的磨耗量为

3.3   曲尖轨磨耗分析结果

车辆通过一次，4个车轮均会与同一尖轨截面发生接触，且4个轮对的运动状态有差别，因此，需要考虑4个车轮与钢轨间磨耗量的累加. 前转向架的前、后轮对左侧车轮编号为1、2，右侧车轮编号为3、4；后转向架的前、后轮对左侧车轮编号为5、6，右侧车轮编号为7、8. 采用两种接触模型得到的磨耗量分布如图6所示，纵坐标X为距离尖轨尖端的距离，从尖轨顶宽为5 mm处开始计算，将两种方法得到的磨耗深度除以最大磨耗深度进行归一化处理.

图  6  尖轨磨耗量分布

Figure  6.  Distribution of wear depth on the switch panel

下载: 全尺寸图片 幻灯片

由图6可知：两种算法得到的磨耗深度分布有明显差异，简化算法得到的最大磨耗深度为精确算法结果的0.75倍. 尖轨靠近基本轨侧简化算法得到的磨耗深度较小，因此无法完全显示. 尽管采用的动力学输入量完全相同，但磨耗深度却出现明显差异，说明两种接触几何算法会对磨耗深度分析造成明显差异，为分析其具体原因，计算沿轨道纵向与曲尖轨接触的4个左侧车轮的接触点纵坐标与接触斑面积，结果如图7所示.

图  7  不同车轮与曲尖轨接触时的接触点纵坐标与接触斑面积

Figure  7.  Longitudinal coordinates of the contact point and contact areas between different wheelsets and the curved switch rail

下载: 全尺寸图片 幻灯片

由图7（a）可知：1） 两种算法得到的车轮1和车轮2接触点纵坐标在轨道前部分几乎一致，但在后半段差异较大；车轮5和车轮6的接触点纵坐标大部分一致，个别位置出现明显差异. 这主要是由于车轮1、2动力学响应较大，当车轮与道岔在基本轨接触时，两种算法差异较小；而车轮与道岔在尖轨接触时，由于简化算法忽略了摇头角与尖轨变截面影响，导致在轮对横移量与摇头角较大时出现明显误差. 2） 精确算法得到的车轮1的接触点纵坐标滞后距离明显小于简化算法，这是由于简化算法考虑等截面尖轨，轮对始终与同一尖轨截面接触；精确算法考虑变截面后，滞后距离越大，则尖轨降低值越大，根据接触点轮轨间垂向距离最小，实际的轮轨接触点滞后距离应比简化算法小. 车轮5、6的摇头角与横移量均较小，车轮与尖轨接触时位置大部分位于尖轨轨顶，尖轨轨顶截面变化较小，因此两种算法得到的接触点纵坐标差异较小.

由图7（b）可知：1） 两种算法得到的接触斑面积差异较小. 这主要是由于动力学仿真分析得到的摇头角较小，在8 mrad以内，因此，两种算法得到的接触点横向位置差异较小，接触点附近的轮轨型面差异较小. 2） 车轮1在尖轨部分的接触斑面积较小，主要是由于车轮横移量较大，车轮与尖轨轨距角接触；车轮2、5、6尖轨部分接触斑面积较大，主要是由于车轮与尖轨轨顶接触. 车轮均在轨道纵向位置3000 ~ 4000 mm处出现接触斑面积骤减，主要是由于此时车轮与尖轨尖端处接触.

为具体分析两种算法得到的磨耗深度差异，计算距离尖轨始端距离为5000 mm处的磨耗量，将磨耗深度除以最大磨耗深度进行无量纲化，其磨耗分布如图8所示.

图  8  钢轨截面磨耗深度分布

Figure  8.  Wear depth on the switch rail cross-section

下载: 全尺寸图片 幻灯片

由图8可知：当Y坐标为670 ~ 680 mm时，简化算法比精确算法得到的最大磨耗深度小50%，这一部分的磨耗主要由第一轮对导致的，第一轮对的横移量与摇头角较大，接触应力较大，因此得到的磨耗深度较大；当Y坐标为685 ~ 710 mm时，这一部分磨耗位于尖轨轨顶，两种算法得到的差异可忽略不计； 仅对第一轮对的接触行为进行分析，第一轮对在距尖轨尖端5000 mm时，摇头角为7.4 mrad，横移量为8.2 mm，法向力为71 kN. 采用两种算法得到的接触斑内切应力分布如图9所示.

图  9  两种算法得到的切应力分布

Figure  9.  Distribution of tangential stress solved by the two methods

下载: 全尺寸图片 幻灯片

由图9可知：两种算法得到的最大切应力量值差异不大，而两种算法得到的切应力分布差异较大，精确算法得到的切应力沿X_c方向显著非对称. 为进一步分析两种算法得到的轮轨切应力分布差异，取X_c = −4.8 mm和Y_c = −1.6 mm处的切应力分析，结果如图10所示.

图  10  切应力沿Y_c和X_c方向分布

Figure  10.  Distribution of tangential stress along Y_c-direction and X_c-direction

下载: 全尺寸图片 幻灯片

由图10可知：当X_c = −4.8 mm时，两种算法得到的应力分布差异较大，简化算法得到的最大切应力为753 MPa，精确算法所得结果为782 MPa；当Y_c =−1.6 mm时，简化算法得到的切应力沿X_c方向对称分布，分布范围为 −6.0 ~ 6.0 mm，精确算法所得结果非对称分布，分布范围为 −3.6 ~ 8.0 mm. 由于两种算法切应力分布存在明显差异，从而磨耗量分布也有明显差异，因此，沿X_c方向求和得到的磨耗深度分布差异较大.

4.   结　论

本文提出了一种更精确的三维非对称接触几何算法，可提高岔区轮轨接触几何计算精确性. 以S1002CN车轮与12号曲尖轨为研究对象，并结合SIMPACK车辆-道岔耦合动力学仿真分析模型和 USFD磨耗模型，分析沿轨道纵向曲尖轨的磨耗分布，结果发现：

1） 精确算法与不考虑尖轨变截面与摇头角的简化算法得到的接触点横坐标差异较小；两种算法得到的接触点纵坐标差异较大. 这是由于精确算法考虑尖轨变截面，尖轨沿轨道纵向宽度与高度均发生变化，根据接触点轮轨垂向距离最小，导致接触点并不位于接触迹线上.

2） 当轮对存在摇头角时，实际的接触斑形状与接触应力沿接触斑纵横向均为非对称分布，而简化算法忽略摇头角与尖轨变截面的影响，导致得到的接触斑结果沿接触斑纵向对称分布，误差较大. 在摇头角为10 mrad，横移量为7.5 mm时，精确算法得到的接触斑面积比简化算法得到的结果大9.2%.

3） 简化算法得到的最大磨耗深度为精确算法结果的0.75倍，这是由于尖轨变截面与摇头角对轮轨接触应力分布与接触斑形状影响较大，尤其是轮缘与轨距角接触时，从而导致简化算法分析尖轨磨耗有较大误差. 本文提出的考虑尖轨变截面的三维非对称接触几何算法可更精确地分析岔区接触几何，为岔区损伤研究提供指导.