基于三维弹性体滚动接触理论的轮轨非平面接触算法

许玉德，严道斌，孙小辉

（1. 同济大学道路与交通工程教育部重点实验室，上海201804；2. 同济大学上海市轨道交通结构耐久与系统安全重点实验室，上海201804；3. 比亚迪汽车工业有限公司，广东深圳518118）

列车的支承、转向、加减速等行为都是由轮轨接触区域的轮轨力完成的，轮轨接触产生的应力是导致轮轨滚动接触疲劳和磨耗的重要因素，特别是在道岔区和曲线地段［1-2］。磨耗的发展导致轮轨接触几何关系恶化，进而影响列车运行的安全性和平稳性。因此，准确计算轮轨接触应力，对于了解轮轨作用机理，研究车辆-轨道耦合作用，评估列车运行的安全性和平稳性有着重要意义。

国内外学者对轮轨接触理论及计算模型开展了深入的研究。Hertz 理论［3］被最早运用于轮轨接触应力计算，但该理论仅适用于平面或近似平面上的椭圆形、单点接触。Kalker［4］提出了三维弹性体滚动接触理论，基于此开发的Contact程序能够较好地计算轮轨非椭圆形、多点接触，是迄今为止使用最为广泛的理论，但该理论使用了平面弹性半空间假设，因此不能较好地解决非平面接触问题。Ayasse等［5］提出了基于虚位移原理的“条带法”计算轮轨法向力以及利用Fastsim 程序［6-7］计算切向力，考虑接触区域曲率对自旋的影响，但仅适用于固定曲率的接触问题。Li［8］基于Kalker的三维弹性体滚动接触理论，提出了利用四分之一空间研究轮轨接触几何的方法，开发的Wear程序适用于共形接触的求解，但影响系数的计算存在精度不足的问题。Vollebregt 等［9］利用有限元方法计算了轮轨共形接触的影响系数，解决了弹性半空间和四分之一空间假设的不足，并将该计算方法加入到Contact 程序中，但这使得Contact程序的计算效率有所下降。Zhu等［10］考虑了轮轨廓形局部曲率变化，并提出了一种改进的半赫兹轮轨接触计算模型。目前，有限元方法是计算轮轨接触应力最为准确的方法，但计算效率过低制约了其用于批量计算的可能性［11］。

在道岔区或曲线地段，钢轨磨耗通常较为严重，轮轨接触面曲率变化较大，导致轮轨非平面接触，进而加剧磨耗，最终对钢轨的使用寿命产生严重影响。准确计算轮轨非平面接触应力，对预测磨耗、优化廓形、延长轮轨使用寿命、保证列车平稳通过［12-14］具有重要意义。基于Kalker 的三维弹性体滚动接触理论，结合轮轨非平面接触几何关系，对最小余能方程中的影响系数进行了修正，提出了一种最小余能方程求解算法，实现了轮轨非平面接触的法向应力、切向应力和黏着区-滑动区的计算。

1 理论方法

1.1 三维弹性体滚动接触理论

如图1 所示，在Kalker 的三维弹性体滚动接触理论中，将可能接触区域划分为矩形单元网格，利用弹性半空间Boussinesq-Cerruti 公式计算得到该区域内任一单元与其他所有单元的位移差，以法向余能最小为目标，通过法向接触算法Norm［4］求解实际接触区域内的法向应力分布；以切向余能最小为目标，通过切向接触算法Tang［4］求解实际接触区域内的切向应力分布以及黏着区-滑动区分布。

三维弹性体滚动接触理论中最小余能方程为

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

式中：C1为余能；Ac为可能接触区域；h为法向间隙；μn为法向位移；pn为法向接触力；wτ为切向刚性滑移量；μτ为切向位移；μ′τ为前一时刻切向位移；pτ为切向接触力；g 为滑动摩擦系数；P 为作用在接触面上的竖向集中力。

对于平面接触问题的求解，该理论认为任一单元的位移与接触面上作用力的对应关系为

式中：A 为可能接触区域各单元的影响系数；Ann为在法向方向作用单位力时产生的法向位移，即法向影响系数；Aττ为在切向方向作用单位力时产生的切向位移，即切向影响系数。

1.2 影响系数的修正

在三维弹性体滚动接触理论中，法向影响系数Ann的计算采用的是弹性半空间的Boussinesq 解析解，切向影响系数Aττ的计算则采用的是Cerruti解析解。在轮轨非平面接触问题中，弹性半空间假设不再成立，这就造成了Kalker的Contact程序无法解决轮轨非平面接触问题，因此有必要对非平面接触时的影响系数进行修正。

参考文献［15］的做法，如图2所示，引入作用力F与曲面的法向角α用于影响系数修正，图中符号含义、具体修正方法及步骤与文献［16］完全相同，并在文献［16］的基础上完善了切向力作用时的影响系数修正公式。

图2 影响系数修正Fig.2 Modification of the influence coefficient

根据文献［15-16］的思路，当可能接触区域内单元J 受到单位作用力时，可能接触区域内任一单元I产生的影响系数可表示为

式中：AIbJa为由Boussinesq-Cerruti 公式计算得到的影响系数，表示在单元J 受a 方向的单位作用力时，单元I 在b 方向的影响系数；为修正后的影响系数；a，b=1，2，3。为方便说明，将法向n、横切向s、纵切向t分别用数字1、2、3表示，如即为，后文中的方向均按此方法表示。

1.3 影响系数矩阵的构造及运用

影响系数修正后，将可能接触区域离散为由面积Δx1×Δx2的矩形构成的网格，设横切向网格数量为M，纵切向网格数量为N，共计MN 个矩形网格。将离散的二维矩形网格从1至MN进行编号，在网格J（J=1，2，···，MN）的a（a=1，2，3）方向作用单位力时，记可能接触区域内所有网格的b（b=1，2，3）方向的影响系数

此时，可能接触区域内所有网格的影响系数可构造为MN×MN维的影响系数矩阵，如下所示：

记沿着a（a=1，2，3）方向作用在可能接触区域所有网格内的接触力

则所有网格在b（b=1，2，3）方向产生的位移可表示为

根据以上论述可知，通过影响系数矩阵构造，对于任意轮轨接触力的作用，利用式（3）～（7）就能计算离散后接触面上所有单元的位移，进一步可以计算该区域内任一单元与其他所有单元的位移差。结合目标函数，可以实现轮轨接触区域法向应力、切向应力和黏着区-滑动区的计算。

2 算法实现

2.1 最小余能方程离散化

在计算轮轨非平面接触时，不再以法向余能最小和切向余能最小为目标函数，而是将总余能最小作为目标函数，求解Kalker的最小余能方程。

将修正后的影响系数代入式（1）中，写成离散形式并展开，所得到的即为轮轨非平面接触条件下的最小余能方程，如下所示：

式中：hI为网格I（I=1，2，…，MN）处法向间隙为网格I 处a（a=1，2，3）方向位移；为网格I 处a（a=1，2，3）方向接触力；为网格I处b（b=2，3）方向刚性滑移量为网格I 处b（b=2，3）方向前一时刻位移；αI为网格I处法向角。

此时，将所有变量表征为向量或矩阵的形式，则在轮轨非平面接触条件下，最小余能方程的离散形式如 下所示：

式 中 ： h 为 法 向 间 隙 矩 阵 ， 即h=[h1h2… hMN]；w 为刚性滑移量矩阵，即w=[w1w2… wMN]。

2.2 最小余能方程求解分析

在式（9）中，p(1)、p(2)、p(3)3组接触力需要求解，从目标函数和约束条件来看，这是一个非线性约束规划问题。首先，将式（9）中约束条件定义如下：

式中：◦为自定义符号，表示相同维度的2 个矩阵之间相同位置元素对应相乘，得到的是相同维度的新矩阵；cos α和sin α分别为法向角余弦值和正弦值矩阵 ， 即 cos α=[cos α1cos α2… cos αMN]，sin α=[sin α1sin α2… sin αMN]。

对于式（10），按照最优解的Kuhn-Tucker 必要条件［17］，必然存在互补变量乘子≥0，≥0，∈R，使得下式同时成立：

式中：∇为偏导符号；λ(a)为互补变量乘子矩阵，即

将式（9）中的C3分别对p(a)I求偏导，即∇C3=，此时式（11）中的L1可全部展开，计算式如下所示：

式（11）、（12）中的6 个方程L11、L12、L13、L2、L3、L4即构成了最优解的Kuhn-Tucker 必要条件展开式，为非线性方程组。6 个方程中变量的数量包括：p(1)I、p(2)I、p(3)I、λ(1)I、λ(2)I各有MN 个，λ(3)I有1 个，总数量为（5MN+1）个。独立方程的数量包括：L11、L12、L13各有MN个，L2、L3各有MN个，L4有1个，总数量为（5MN+1）个。变量与独立方程的数量相同，因此可以进行最优化求解。

2.3 最小余能方程求解算法

直接采用牛顿-拉夫逊方法（N-R 法）［18］对6个方程进行最优求解是困难的，因为方程L2是一个三次方程，而N-R 法只是二次收敛的。对于类似的=0 互补松弛方程，和中必然有一个因子为0，因此可以通过判别，使其中一个因子为0，从而将方程降为二次方程，同时将方程简化。

将整个计算区域分为接触区域和非接触区域。对于计算区域内的任一网格I（I=1，2，…，MN），在接触区域内的所有法向接触力必然有＞0，同时又要满足互补松弛方程L3=0，因此在接触区域内对应的=0。在非接触区域内，又有=0，则＞0。因此，当假定了接触区域后，就可以将互补松弛方程L3释放掉。

对于接触区域，又分为黏着区和滑动区。对于在接触区域内的任一网格I，在黏着区内必然有，同时又要满足互补松 弛 方 程，因此在黏着区内对应的=0。在滑动区内，又有，则。通过 以上判定将L2降为二次方程，就能较好地使用N-R法。算法的基本步骤如下所示：

步骤1对可能接触区域进行网格划分，输入已知参数，以此计算影响系数、法向间隙等基本参数。

步骤2假定可能接触区域全部为接触区域，并设置为黏着区，此时的初始解均为0。

步骤3计算式（11）中的线性方程L1和L4，得到的初始解。

步骤4检查接触区域＞0是否成立。若是，则执行步骤（5）；若不是，则将≤0的单元设定为非接触区域，执行步骤（8）。

步骤5检查非接触区域＞0 是否成立。若是，则执行步骤6；若不是，则将≤0 的单元设定为接触区域，执行步骤8。

步 骤 6检 查 黏 着 区是否成立。若是，则执行步骤7；若不是，则将的单元设定为滑动区，同时设定切向接触力

步骤7检查滑动区＞0是否成立。若是，则当前解即为计算结果，计算结束；若不是，则将≤0的单元设定为黏着区，执行步骤8。

步骤8以上步骤中的作为初值T0，首先计算Tk（k=0，1，2，…）对应的函数值L(Tk)和雅克比矩阵JL(Tk)，其次按照迭代公式ΔTk=-JL(Tk)-1L(Tk)计算迭代步长。若|ΔTk|≤ξ成立，其中ξ为精度条件，则执行步骤3；若不成立，则更新变量使Tk+1=Tk+ΔTk，重新执行步骤8，直到|ΔTk|≤ξ成立为止。

3 仿真验证

为验证所提算法的准确性，以LM 型车轮和CHN75钢轨为对象，利用所提算法计算轮轨轨距角非平面接触时的法向应力、切向应力和黏着区-滑动区。同时，在有限元软件中建立完全一致的工况并计算对应的结果，通过对2 种算法计算结果的比较来判断所提算法的准确性。设置的工况及对应的参数如表1 所示。需要说明的是，不考虑钢轨截面纵向的变化。

表1 验证参数Tab.1 Parameters used for validation

如图3 所示，在横移量为+5.5 mm 条件下，车轮与钢轨发生轨距角非平面接触，该接触为两点接触，计算结果可用于验证。

图3 轮轨接触二维断面Fig.3 Two-dimensional cross section of wheel-rail contact

3.1 法向应力

如图4所示，从法向应力计算结果来看，所提算法与有限元方法得到的结果是接近的。所提算法相比于有限元方法，在“高峰”处最大值差异为-2.6%，在“低峰”处最大值差异为-6.0%。

3.2 切向应力

图4 法向应力结果Fig.4 Results of normal stress

如图5所示，从切向应力计算结果来看，所提算法与有限元方法得到的结果也是接近的。所提算法相比于有限元方法，在“高峰”处最大值差异为-2.0%，在“低峰”处最大值差异为-6.0%。

图5 切向应力结果Fig.5 Results of shear stress

3.3 黏着区-滑动区

如图6 所示，从黏着区-滑动区分布的情况来看，所提算法相比于有限元方法，在“高峰”处的分布是比较接近的。值得注意的是，在“低峰”处有限元方法计算结果中出现了黏着区，但所提算法计算结果中并未出现。进一步对黏着区-滑动区的面积进行了对比，结果如表2 所示。所提算法相比于有限元方法，差异基本都控制在了7%以内，这种差异还是可以接受的。

图6 黏着区-滑动区分布Fig.6 Distribution of stick-slip region

表2 黏着区-滑动区面积Tab.2 Area of stick-slip region 单位：mm2

综合上述对比，可以认为所提算法能够有效计算轮轨非平面接触时的法向应力、切向应力和黏着区-滑动区。需要指出的是，由于有限元方法计算耗时过长，在计算本验证工况时超过1 h，而所提算法计算时长极短（不超过1 s），在计算效率上，所提算法更有优势，对车辆-轨道耦合动力学这类需要批量计算轮轨非平面接触问题的求解也更为有利。

4 算例

利用所提算法研究钢轨在不同磨耗状态下的轮轨接触特性。以我国某重载铁路开行的27 t轴重列车为例进行计算，选取通过总重约为60、90、200 MGT（1 MGT=106t）的磨耗钢轨，取自同一曲线外轨测点。磨耗钢轨型面如图7所示，计算参数如表3所示。

图7 磨耗钢轨型面Fig.7 Worn rail profiles

表3 算例参数Tab.3 Parameters used for calculation

需要说明的是，第2、3次测量间隔时间较长，因此本算例仅是对不同磨耗状态下的轮轨接触特性变化趋势做出判断和分析。

4.1 法向应力

如图8 所示，在+13.0 mm 的横移量条件下，随着磨耗的增加，在60 MGT 时为典型的两点接触。虽然在90 MGT 时仍然属于两点接触，但是由明显的“双峰”曲面向“单峰”曲面变化。到200 MGT时，接触状态过渡至共形接触，应力峰值出现在接触中心。在此过程中，接触斑形状从60 MGT 时的纵向双“椭圆形”，转变至90 MGT时的横向“葫芦形”，直至200 MGT 时变成狭长的横向“椭圆形”。在此过程中，最大法向应力在逐渐减小，而接触面积则逐渐增大。

图8 法向应力演变Fig.8 Evolution of normal stress

4.2 切向应力

如图9 所示，在+13.0 mm 横移量的条件下，最大切向应力在逐渐减小，值得关注的是，切向应力峰值点的位置由60 MGT和90 MGT的双接触中心逐渐演变至200 MGT的多接触中心。

图9 切向应力演变Fig.9 Evolution of shear stress

4.3 黏着区-滑动区

如图10 所示，在+13.0 mm 横移量条件下，黏着区-滑动区分布形状出现了较为显著的变化，黏着区的形状由纵向“椭圆形”逐渐过渡至横向“椭圆形”。表4列出了3种磨耗状态下黏着区和滑动区的面积。可以看出，随着总面积的增大，黏着区和滑动区的面积也逐渐增大。

图10 黏着区-滑动区演变Fig.10 Evolution of stick-slip region

表4 算例的黏着区-滑动区面积Tab.4 Area of stick-slip region for calculation

5 结语

基于Kalker 的三维弹性体滚动接触理论，结合轮轨非平面接触几何关系，提出了轮轨非平面接触状态下的影响系数计算公式，并构造了影响系数矩阵，用于最小余能方程的离散化。将离散化最小余能方程的求解转化为非线性规划问题，根据对应非线性方程组最优解的Kuhn-Tucker 必要条件，提出了利用牛顿-拉夫逊方法求最小余能方程最优解的算法，实现了轮轨非平面接触状态下法向应力、切向应力和黏着区-滑动区的计算。通过有限元方法验证了所提算法的准确性，并利用该算法初步探索了钢轨在不同磨耗状态下的轮轨接触特性。

该研究成果为精确、快速计算轮轨非平面接触特性提供了新思路，未来将进一步探索在轮轨磨耗预测、轮轨廓形优化、轮轨使用寿命延长、列车运行安全性和平稳性评估等方面的适用性。