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

许玉德，严道斌，孙小辉; XU　Yude; YAN　Daobin; SUN　Xiaohui

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

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

许玉德 ^1,2
严道斌 ^1,2
孙小辉 ³

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

中图分类号： U211.5

最近更新：2020-03-27

DOI：10.11908/j.issn.0253-374x.19314

摘要

基于Kalker的三维弹性体滚动接触理论，结合轮轨非平面接触几何关系，提出了最小余能方程中影响系数的修正公式。考虑轮轨非平面接触时法向与切向存在的相互作用，对最小余能方程进行离散化。以总余能最小为目标，将离散方程的求解转化为非线性规划问题，并提出了求解算法。利用有限元方法验证了所提算法的准确性。最后，研究了钢轨在不同磨耗状态下的轮轨接触特性。结果表明，在不同磨耗、相同横移量条件下，随着磨耗的增加，轮轨的轨距角接触由两点接触过渡至共形接触，最大接触应力减小，接触斑变得狭长，接触面积增大。

关键词

轮轨接触; 非平面接触; 最小余能方程

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

国内外学者对轮轨接触理论及计算模型开展了深入的研究。Hertz理论^[

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

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

1 理论方法

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

如图1所示，在Kalker的三维弹性体滚动接触理论中，将可能接触区域划分为矩形单元网格，利用弹性半空间Boussinesq‒Cerruti公式计算得到该区域内任一单元与其他所有单元的位移差，以法向余能最小为目标，通过法向接触算法Norm^[

4]求解实际接触区域内的法向应力分布；以切向余能最小为目标，通过切向接触算法Tang^{[参考文献 4
百度学术}4]求解实际接触区域内的切向应力分布以及黏着区‒滑动区分布。

图 1 可能接触区域网格划分

Fig. 1 Meshing of the possible contact region

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

$\begin{array}{l} m i n C_{1} = \int_{A_{c}} (h + \frac{1}{2} μ_{n}) p_{n} d S + \\ \int_{A_{c}} (w_{τ} + \frac{1}{2} μ_{τ} - μ_{τ}^{'}) p_{τ} d S \end{array}$

$s . t . g p_{n} - |p_{τ}| \geq 0$ ，

$p_{n} \geq 0$ ,

$\sum p_{n} = P$

(1)

式中： $C$ ₁为余能； $A_{c}$ 为可能接触区域； $h$ 为法向间隙； $μ_{n}$ 为法向位移； $p_{n}$ 为法向接触力； $w_{τ}$ 为切向刚性滑移量； $μ_{τ}$ 为切向位移； $μ_{τ}^{'}$ 为前一时刻切向位移； $p_{τ}$ 为切向接触力； $g$ 为滑动摩擦系数； $P$ 为作用在接触面上的竖向集中力。

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

$μ_{n} = p_{n} A_{n n} ， μ_{τ} = p_{τ} A_{τ τ}$

(2)

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

1.2　影响系数的修正

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

参考文献[

15]的做法，如图2所示，引入作用力

$F$ 与曲面的法向角

$α$ 用于影响系数修正，图中符号含义、具体修正方法及步骤与文献[16]完全相同，并在文献[16]的基础上完善了切向力作用时的影响系数修正公式。

图 2 影响系数修正

Fig. 2 Modification of the influence coefficient

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

$\{\begin{array}{l} A_{I, J}^{(b 1)} = A_{I b J 1} c o s α - A_{I b J 2} s i n α \\ A_{I, J}^{(b 2)} = A_{I b J 2} c o s α + A_{I b J 1} s i n α \\ A_{I, J}^{(b 3)} = A_{I b J 3} \end{array}$

(3)

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

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

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

$d_{J}^{(b a)} = [\begin{matrix} A_{1, J}^{(b a)} & A_{2, J}^{(b a)} & \dots & A_{M N, J}^{(b a)} \end{matrix}]$

(4)

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

$A_{}^{(b a)} = [\begin{matrix} d_{1}^{(b a)} \\ d_{2}^{(b a)} \\ ⋮ \\ d_{M N}^{(b a)} \end{matrix}] = [\begin{matrix} A_{1,1}^{(b a)} & A_{2,1}^{(b a)} & \dots & A_{M N, 1}^{(b a)} \\ A_{1,2}^{(b a)} & A_{2,2}^{(b a)} & \dots & A_{M N, 2}^{(b a)} \\ ⋮ & ⋮ & ⋮ \\ A_{1, M N}^{(b a)} & A_{2, M N}^{(b a)} & \dots & A_{M N, M N}^{(b a)} \end{matrix}]$

(5)

$p_{}^{(a)} = [\begin{matrix} p_{1}^{(a)} & p_{2}^{(a)} & \dots & p_{M N}^{(a)} \end{matrix}]$

(6)

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

$μ_{}^{(b)} = [μ_{1}^{(b)} μ_{2}^{(b)} \dots μ_{M N}^{(b)}] = \overset{3}{\sum_{a = 1}} p_{}^{(a)} A^{(b a)}$

(7)

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

2 算法实现

2.1　最小余能方程离散化

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

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

$m i n C_{2} = \overset{M N}{\sum_{I = 1}} ((h_{I} + \frac{1}{2} μ_{I}^{(1)}) p_{I}^{(1)} + (w_{I}^{(2)} + \frac{1}{2} μ_{I}^{(2)} - {μ_{I}^{'}}^{(2)}) p_{I}^{(2)} + (w_{I}^{(3)} + \frac{1}{2} μ_{I}^{(3)} - {μ_{I}^{'}}^{(3)}) p_{I}^{(3)})$

$s . t . {(g p_{I}^{(1)})}^{2} - ({(p_{I}^{(2)})}^{2} + {(p_{I}^{(3)})}^{2}) \geq 0$ ,

$p_{I}^{(1)} \geq 0$ ，

$\overset{M N}{\sum_{I = 1}} (p_{I}^{(1)} c o s α_{I} + p_{I}^{(2)} s i n α_{I}) = P$

(8)

式中： $h_{I}$ 为网格I（ $I = 1 ， 2 ， \dots ， M N$ ）处法向间隙； $μ_{I}^{(a)}$ 为网格I处 $a$ （ $a = 1 ， 2 ， 3$ ）方向位移； $p_{I}^{(a)}$ 为网格I处 $a$ （ $a = 1 ， 2 ， 3$ ）方向接触力； $w_{I}^{(b)}$ 为网格I处 $b$ （ $b = 2 ， 3$ ）方向刚性滑移量； ${μ_{I}^{'}}^{(b)}$ 为网格I处 $b$ （ $b = 2 ， 3$ ）方向前一时刻位移； $α_{I}$ 为网格I处法向角。

对于前一时刻位移 ${μ_{I}^{'}}^{(b)}$ ，仅考虑车轮稳态滚动。在稳态滚动的过程中，轮轨接触力保持不变，前一时刻接触力 $p_{I}^{'}$ 与当前时刻接触力 $p_{I}$ 的大小和分布一致，仅是在离散的网格中向前整体移动了一个网格的距离。因此，在计算影响系数时，需要对可能接触区域进行拓宽，从而将接触前沿和接触后沿包括进来，此时可能接触区域网格数增至 $M (N + 1)$ ，记拓宽后的影响系数矩阵为 $\overset{̿}{A}$ ，其维度为 $M (N + 1) \times M (N + 1)$ 。记接触前沿网格编号的集合为X，在前一时刻，接触前沿在可能接触区域外，记 $\bar{A} = \overset{̿}{A}| X$ 。此时，将所有变量表征为向量或矩阵的形式，则在轮轨非平面接触条件下，最小余能方程的离散形式如下所示：

$m i n C_{3} = h p^{(1) T} + \overset{3}{\sum_{a = 1}} \frac{1}{2} p^{(a)} A^{(1 a)} p^{(1) T} + \overset{3}{\sum_{b = 2}} w^{(b)} p^{(b) T} + \overset{3}{\sum_{a = 1}} \overset{3}{\sum_{b = 2}} (\frac{1}{2} p^{(a)} (A^{(b a)} - 2 {\bar{A}}^{(b a)}) p^{(b) T})$

s.t.

${(g p_{I}^{(1)})}^{2} - ({(p_{I}^{(2)})}^{2} + {(p_{I}^{(3)})}^{2}) \geq 0,$

$p_{I}^{(1)} \geq 0$ ,

$\overset{M N}{\sum_{I = 1}} (p_{I}^{(1)} c o s α_{I} + p_{I}^{(2)} s i n α_{I}) = P$

(9)

式中： $h$ 为法向间隙矩阵，即 $h = [\begin{matrix} h_{1} & h_{2} & \dots & h_{M N} \end{matrix}]$ ； $w$ 为刚性滑移量矩阵，即 $w = [\begin{matrix} w_{1} & w_{2} & \dots & w_{M N} \end{matrix}]$ 。

2.2　最小余能方程求解分析

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

$\{\begin{array}{l} G_{1} = g^{2} (p_{}^{(1)} \circ p_{}^{(1)}) - ((p_{}^{(2)} \circ p_{}^{(2)}) + (p_{}^{(3)} \circ p_{}^{(3)})) \\ G_{2} = p^{(1)} \\ G_{3} = p^{(1)} \circ (c o s α) + p^{(2)} \circ (s i n α) \end{array}$

(10)

式中： $\circ$ 为自定义符号，表示相同维度的2个矩阵之间相同位置元素对应相乘，得到的是相同维度的新矩阵； $c o s α$ 和 $s i n α$ 分别为法向角余弦值和正弦值矩阵，即 $c o s α = [\begin{matrix} c o s α_{1} & c o s α_{2} & \dots & c o s α_{M N} \end{matrix}]$ ， $s i n α = [\begin{matrix} s i n α_{1} & s i n α_{2} & \dots & s i n α_{M N} \end{matrix}]$ 。

对于式(10)，按照最优解的Kuhn‒Tucker必要条件^[

17]，必然存在互补变量乘子

$λ_{I}^{(1)} \geq 0$ ，

$λ_{I}^{(2)} \geq 0$ ，

$λ_{I}^{(3)} \in$ R，使得下式同时成立：

$\{\begin{array}{l} L_{1} : \nabla C_{3} - λ^{(1)} \circ (\nabla G_{1}) - λ^{(2)} \circ (\nabla G_{2}) - \\ λ^{(3)} \circ (\nabla G_{3}) = 0 \\ L_{2} : λ_{I}^{(1)} ({(g p_{I}^{(1)})}^{2} - ({(p_{I}^{(2)})}^{2} + {(p_{I}^{(3)})}^{2})) = 0 \\ L_{3} : λ_{I}^{(2)} p_{I}^{(1)} = 0 \\ L_{4} : λ_{I}^{(3)} P = 0 \end{array}$

(11)

式中： $\nabla$ 为偏导符号； $λ^{(a)}$ 为互补变量乘子矩阵，即 $λ^{(a)} = [\begin{matrix} λ_{1}^{(a)} & λ_{2}^{(a)} & \dots & λ_{M N}^{(a)} \end{matrix}]$ 。

将式(9)中的 $C_{3}$ 分别对 $p_{I}^{(a)}$ 求偏导，即 $\nabla C_{3} = \frac{\partial C_{3}}{\partial p_{I}^{(a)}}$ ，此时式(11)中的 $L_{1}$ 可全部展开,计算式如下所示：

$\{\begin{array}{l} L_{11} : h + \frac{1}{2} p^{(1)} (A^{(11)} + A^{(11) T}) + \frac{1}{2} p^{(2)} (A^{(12)} + A^{(21) T} - 2 {\bar{A}}^{(21) T}) + \frac{1}{2} p^{(3)} (A^{(13)} + A^{(31) T} - 2 {\bar{A}}^{(31) T}) - \\ 2 g^{2} λ^{(1)} \circ p^{(1)} - λ^{(2)} - λ^{(3)} \circ (c o s α) = 0 \\ L_{12} : w^{(2)} + \frac{1}{2} p^{(1)} (A^{(21)} - 2 {\bar{A}}^{(21)} + A^{(12) T}) + \frac{1}{2} p^{(2)} (A^{(22)} - 2 {\bar{A}}^{(22)} + A^{(22) T} - {\bar{A}}^{(22) T}) + \\ \frac{1}{2} p^{(3)} (A^{(23)} - 2 {\bar{A}}^{(23)} + A^{(32) T} - 2 {\bar{A}}^{(32) T}) + 2 λ^{(1)} \circ p^{(2)} - λ^{(3)} \circ (s i n α) = 0 \\ L_{13} : w^{(3)} + \frac{1}{2} p^{(1)} (A^{(31)} - 2 {\bar{A}}^{(31)} + A^{(13) T}) + \frac{1}{2} p^{(2)} (A^{(32)} - 2 {\bar{A}}^{(32)} + A^{(23) T} - 2 {\bar{A}}^{(23) T}) + \\ \frac{1}{2} p^{(3)} (A^{(33)} - 2 {\bar{A}}^{(33)} + A^{(33) T} - 2 {\bar{A}}^{(33) T}) + 2 λ^{(1)} \circ p^{(3)} = 0 \end{array}$

(12)

式(11)、(12)中的6个方程 $L_{11}$ 、 $L_{12}$ 、 $L_{13}$ 、 $L_{2}$ 、 $L_{3}$ 、 $L_{4}$ 即构成了最优解的Kuhn‒Tucker必要条件展开式，为非线性方程组。6个方程中变量的数量包括： $p_{I}^{(1)} 、 p_{I}^{(2)} 、 p_{I}^{(3)} 、 λ_{I}^{(1)} 、 λ_{I}^{(2)}$ 各有 $M N$ 个， $λ_{I}^{(3)}$ 有1个，总数量为（ $5 M N + 1$ ）个。独立方程的数量包括： $L_{11}$ 、 $L_{12}$ 、 $L_{13}$ 各有 $M N$ 个， $L_{2}$ 、 $L_{3}$ 各有 $M N$ 个， $L_{4}$ 有1个，总数量为（ $5 M N + 1$ ）个。变量与独立方程的数量相同，因此可以进行最优化求解。

2.3　最小余能方程求解算法

直接采用牛顿‒拉夫逊方法（N‒R法）^[

18]对6个方程进行最优求解是困难的，因为方程

$L_{2}$ 是一个三次方程，而N‒R法只是二次收敛的。对于类似的

$λ_{I}^{(2)} p_{I}^{(1)} = 0$ 互补松弛方程，

$λ_{I}^{(2)}$ 和

$p_{I}^{(1)}$ 中必然有一个因子为0，因此可以通过判别，使其中一个因子为0，从而将方程降为二次方程，同时将方程简化。

将整个计算区域分为接触区域和非接触区域。对于计算区域内的任一网格I（ $I = 1,2 ， \dots ， M N$ ），在接触区域内的所有法向接触力必然有 $p_{I}^{(1)} > 0$ ，同时又要满足互补松弛方程 $L_{3} = 0$ ，因此在接触区域内对应的 $λ_{I}^{(2)} = 0$ 。在非接触区域内，又有 $p_{I}^{(1)} = 0$ ，则 $λ_{I}^{(2)} > 0$ 。因此，当假定了接触区域后，就可以将互补松弛方程 $L_{3}$ 释放掉。

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

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

步骤2 假定可能接触区域全部为接触区域，并设置为黏着区，此时 $λ_{I}^{(1)}$ 、 $λ_{I}^{(2)}$ 的初始解均为0。

步骤3 计算式(11)中的线性方程 $L_{1}$ 和 $L_{4}$ ，得到 $p_{I}^{(1)}$ 、 $p_{I}^{(2)}$ 、 $p_{I}^{(3)}$ 、 $λ_{I}^{(3)}$ 的初始解。

步骤4 检查接触区域 $p_{I}^{(1)} > 0$ 是否成立。若是，则执行步骤（5）；若不是，则将 $p_{I}^{(1)} \leq 0$ 的单元设定为非接触区域，执行步骤（8）。

步骤5 检查非接触区域 $λ_{I}^{(2)} > 0$ 是否成立。若是，则执行步骤6；若不是，则将 $λ_{I}^{(2)} \leq 0$ 的单元设定为接触区域，执行步骤8。

步骤6 检查黏着区 ${(g p_{I}^{(1)})}^{2} - ({(p_{I}^{(2)})}^{2} +$

${(p_{I}^{(3)})}^{2}) > 0$ 是否成立。若是，则执行步骤7；若不是，则将 ${(g p_{I}^{(1)})}^{2} - ({(p_{I}^{(2)})}^{2} + {(p_{I}^{(3)})}^{2}) \leq 0$ 的单元设定为滑动区，同时设定切向接触力

$p_{I}^{(a)} = g p_{I}^{(1)} p_{I}^{(a)} {({(p_{I}^{(2)})}^{2} + {(p_{I}^{(3)})}^{2})}^{- 1} ， a = 2 ， 3$

(13)

步骤7 检查滑动区 $λ_{I}^{(1)} > 0$ 是否成立。若是，则当前解即为计算结果，计算结束；若不是，则将 $λ_{I}^{(1)} \leq 0$ 的单元设定为黏着区，执行步骤8。

步骤8 以上步骤中的 $p_{I}^{(1)}$ 、 $p_{I}^{(2)}$ 、 $p_{I}^{(3)}$ 、 $λ_{I}^{(1)}$ 、 $λ_{I}^{(2)}$ 、 $λ_{I}^{(3)}$ 作为初值 $T_{0}$ ，首先计算 $T_{k}$ （ $k = 0 ， 1 ， 2, \dots$ ）对应的函数值 $L (T_{k})$ 和雅克比矩阵 $J_{L} (T_{k})$ ，其次按照迭代公式 $Δ T_{k} = - J_{L} (T_{k})^{- 1} L (T_{k})$ 计算迭代步长。若 $|Δ T_{k}| \leq ξ$ 成立，其中 $ξ$ 为精度条件，则执行步骤3；若不成立，则更新变量使 $T_{k + 1} = T_{k} + Δ T_{k}$ ，重新执行步骤8，直到 $|Δ T_{k}| \leq ξ$ 成立为止。

3 仿真验证

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

表 1 验证参数

Tab.1 Parameters used for validation

参数	数值
车轮型面	LM
车辆轴重/t	30
车辆速度/（km·h^-1）	60
钢轨型面	CHN75
弹性模量/GPa	210
泊松比	0.28
轨底坡	1∶40
滑动摩擦系数	0.3
横移量/mm	+5.5
网格尺寸/（mm×mm）	1.0×1.0

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

图3 轮轨接触二维断面

Fig.3 Two-dimensional cross section of wheel-rail contact

3.1　法向应力

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

图4 法向应力结果

Fig.4 Results of normal stress

3.2　切向应力

如图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 ( mm² )

方法	“高峰”		“低峰”		合计面积
方法	黏着区面积	滑动区面积	黏着区面积	滑动区面积	合计面积
有限元方法	86	80	5	30	201
所提算法	81	76	0	31	187
差异/%	-5.8	-5.0		+3.3	-7.0

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

4 算例

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

图7 磨耗钢轨型面

Fig.7 Worn rail profiles

表3 算例参数

Tab.3 Parameters used for calculation

参数	数值
车轮型面	LM
车辆轴重/t	27
车辆速度/（km·h^-1）	60
弹性模量/GPa	210
泊松比	0.28
轨底坡	1∶40
滑动摩擦系数	0.3
横移量/mm	+13.0
网格尺寸/（mm×mm）	1.0×1.0

需要说明的是，第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

通过总重/MGT	黏着区面积/mm²	滑动区面积/mm²	合计面积/mm²
60	151	128	279
90	215	169	384
200	263	234	497

5 结语

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

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

参考文献

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. [百度学术]

MEYMAND S, KEYLIN A, AHMADIAN M. A survey of wheel-rail contact models for rail vehicles[J]. Vehicle System Dynamics, 2016, 54:386. [百度学术]

HERTZ H. Über die berührung fester elastischer körper[J]. Reine und Angewandte Mathematik, 1882, 92:156. [百度学术]

KALKER J. Three-dimensional elastic bodies in rolling contact[M]. Dordrecht: Kluwer Academic Publishers, 1990. [百度学术]

AYASSE J, CHOLLET H. Determination of the wheel rail contact patch in semi-Hertzian conditions[J]. Vehicle System Dynamics, 2005, 43:161. [百度学术]

KALKER J. A fast algorithm for the simplified theory of rolling contact[J]. Vehicle System Dynamics, 1982, 11:1. [百度学术]

KALKER J. A simplified theory for non-Hertzian contact[J]. Vehicle System Dynamics, 1983, 12:43. [百度学术]

LI Z. Wheel-rail rolling contact and its application to wear simulation[D]. Delft: Delft University of Technology, 2002. [百度学术]

VOLLEBREGT E, SEGAL G. Solving conformal wheel-rail rolling contact problems[J]. Vehicle System Dynamics, 2014, 52(S1):455. [百度学术]

ZHU B, ZENG J, ZHANG D, et al. A non-Hertzian wheel-rail contact model considering wheelset yaw and its application in wheel wear prediction[J]. Wear, 2019, 432:1. [百度学术]

MA Y, MARKINE V, MASHAL A, et al. Modelling verification and influence of operational patterns on tribological behaviour of wheel-rail interaction[J]. Tribology International, 2017, 114:264. [百度学术]

ZHOU Y, HAN Y, MU D, et al. Prediction of the coexistence of rail head check initiation and wear growth[J]. International Journal of Fatigue, 2018, 112: 289. [百度学术]

BUTINI E, MARINI L, MEACCI M, et al. An innovative model for the prediction of wheel-rail wear and rolling contact fatigue[J]. Wear, 2019, 436:1. [百度学术]

FIRLIK B, STASKIEWICZ T, JASKOWSKI W, et al. Optimization of a tram wheel profile using a biologically inspired algorithm[J]. Wear, 2019, 430:12. [百度学术]

BLANCO L J, SANTAMARIA J, VADILLO E, et al. On the influence of conformity on wheel-rail rolling contact mechanics[J]. Tribology International, 2016, 103:647. [百度学术]

许玉德, 严道斌, 孙小辉, 等. 重载铁路钢轨磨耗状态下的轮轨法向接触特性[J]. 同济大学学报：自然科学版, 2019, 47(5):663. [百度学术]

XU Yude, YAN Daobin, SUN Xiaohui, et al. Wheel-rail normal contact characteristics of heavy haul worn rails[J]. Journal of Tongji University：Natural Science, 2019, 47(5):663. [百度学术]

李向利, 刘红卫, 黄亚魁. 求解特定线性互补问题的牛顿KKT内点法[J]. 应用数学学报, 2010, 33(5):889. [百度学术]

LI Xiangli, LIU Hongwei, HUANG Yakui. Newton-KKT interior-point methods for special linear complementarity problem[J]. Acta Mathematicae Applicatae Sinica, 2010, 33(5):889. [百度学术]

冯冬冬. 一类精细修正牛顿法和拟牛顿法研究[D].长沙: 中南大学, 2011. [百度学术]

FENG Dongdong. Research on nonmonotone modified Newton and quasi-Newton methods[D]. Changsha: Central South University, 2011. [百度学术]

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

摘要

关键词

1 理论方法

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

1.2 影响系数的修正

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

2 算法实现

2.1 最小余能方程离散化

2.2 最小余能方程求解分析

2.3 最小余能方程求解算法