崔 潇，姚建伟，孙丽霞

(1.中国铁道科学研究院 研究生部，北京 100081；2.中国铁道科学研究院集团有限公司 铁道科学技术研究发展中心，北京 100081)

列车在高速运行过程中，轮轨接触非稳态特性逐渐明显，钢轨波浪形磨耗、车轮多边形磨损、轮轨接触噪声等都与其密切相关。轮轨非稳态接触是指当车轮与钢轨发生滚动接触时，接触表面上的任意质点从进入接触斑前沿到离开接触斑后沿的过程中，轮轨接触表面几何外形、轮轨间法向力、蠕滑率等接触相关参数发生较快变化的现象。Kalker最先研究了轮轨二维非稳态问题[1-2]，利用复变函数方法将非稳态问题转换为实轴上的Riemann-Hilbert问题，并通过数值积分进行求解。Gonzales等[3]采用边界元方法计算了二维非稳态问题，获得了与Kalker一致的结果。此后，学者们将研究范围拓展到三维情况。Kalker基于边界元方法和最小余能原理开发了CONTACT程序，并由Vollebregt[4]及其他学者逐步扩展。CONTACT程序目前已经成为轮轨非稳态接触计算中公认的参考标准。CONTACT程序虽然计算结果精确，但计算量较大，无法在动力学软件中直接集成使用。因此，很多学者进行了非稳态接触简化模型的研究。

1 计算理论

1.1 轮轨非稳态滚动接触计算理论

利用Kalker三维滚动接触精确解方法，计算轮轨非稳态滚动接触状态。基于最小余能原理的CONTACT程序的计算理论可以表示为

(1)

最小值的求解可以视为一个优化问题。当2个接触体选用的材料参数不同，无法视为准同一问题时，该问题的求解可以分为法向求解和切向求解分别迭代计算。当可以视为准同一的情况下，则可以使用KOMBI算法直接迭代求解。

1.2 轮轨非稳态滚动接触蠕滑率/蠕滑力传递函数

轮轨非稳态滚动接触是一个复杂的问题，它与法向力、蠕滑率、接触几何、摩擦因数等诸多因素相关。基于CONTACT程序，将非稳态参数(如蠕滑率)设为简谐波动形式，计算蠕滑力相位与幅值增益。然后对非稳态滚动接触系统进行传递函数估计，从而评估非稳态滚动接触模型单个影响因素的性质。

设系统的传递频率响应函数G(jω)为

(2)

式中：a0,a1,…,an为传递函数的系数；ω为蠕滑率简谐波动的频率。

利用倒幅相特性并将其分别展开为实部与虚部：

(3)

(4)

(5)

将计算获得的频率、幅值增益与相位代入到式(4)、式(5)，再通过最小二乘法，获得式(4)、式(5)的系数，最终确定传递函数。

2 计算工况及结果分析

选取我国高速铁路LMA型踏面车轮和CHN60型钢轨。车轮半径为430 mm，接触点处车轮主轮廓线半径为450 mm，接触点处钢轨头主轮廓线半径为300 mm。计算可知，接触斑长轴半径a≈6.61 mm，短轴半径b≈3.71 mm，则静态接触斑的长短轴比a/b≈1.78。动静摩擦因数均取0.4,静态法向力为70 kN。

在仅考虑单一方向简谐波动蠕滑率输入的工况下，输入的蠕滑率ξ采用如下形式

(6)

式中：x为前进方向距离；ξ0为参考蠕滑率；Δξ为蠕滑率波动分量；Lξ为蠕滑率波长；φ0为初始相位角。

根据Kalker线性理论，Δξ和蠕滑力波动分量ΔFξ之间的幅值增益gξ为

(7)

式中：G为剪切模量；C为Kalker线性理论常数。

Δξ和ΔFξ之间的相位差为

Δφ=φF-φ0

(8)

式中：φF为蠕滑力的相位角。

非稳态蠕滑力相对蠕滑率的滞后效应见图1。

图1 蠕滑力滞后效应

在考虑其他方向蠕滑率共同影响的情况下，使用简谐波动蠕滑率输入时，其他方向蠕滑率采用恒定参考蠕滑率。对同一蠕滑率工况，通过计算多个不同波长比Lξ/a，可以获得一系列的幅值增益和相位。根据公式(9)—式(11)，可将一系列幅值增益和相位的观测值点对(gξ，Δφ)在复平面上采用Nyquist图表示出来。

G(jω)=Re[G(jω)]+j Im[G(jω)]

(9)

A(ω)=|G(jω)|

(10)

(11)

式中：A(ω),φ(ω)分别为传递函数的幅值增益和相位。

采用最小二乘法在复平面上对观测值进行拟合，则可以获得蠕滑率/蠕滑力的传递函数。

2.1 单一方向蠕滑率对传递函数的影响

输入简谐波动的纵向蠕滑率，利用CONTACT程序计算输出纵向蠕滑力。分别选取参考蠕滑率为0.1‰,1.0‰,2.0‰和3.0‰的工况，通过输入多组不同波长的纵向蠕滑率，计算并拟合出纵向蠕滑率与纵向蠕滑力的传递函数。横向蠕滑率/横向蠕滑力、自旋蠕滑率/横向蠕滑力拟合传递函数的方法与其一致。图2为蠕滑力随波长比Lξ/a变化的Nyquist图。

参考蠕滑率/‰ 0.1 1.0 2.0 3.0图2 蠕滑力随波长比Lξ/a变化的Nyquist图

Nyquist图上的横坐标为A(ω)sin[φ(ω)]，即幅值增益在虚轴上的分量;纵坐标为A(ω)cos[φ(ω)]，即幅值增益在实轴上的分量，其中ω取4πa/Lξ。对于图2中同一参考蠕滑率下的每个数据点，输入波长都采用了不同的波长比Lξ/a。图2(a)表示以纵向蠕滑率为输入，纵向蠕滑力作为输出计算幅值增益和相位。

纵向蠕滑率/纵向蠕滑力、横向蠕滑率/横向蠕滑力传递函数可采用一阶传递函数1/[a0+a1(jω)]进行拟合；自旋蠕滑率/横向蠕滑力可采用三阶传递函数1/[a0+a1(jω)+a2(jω)2+a3(jω)3]进行拟合。传递函数系数如表1—表3所示。

表1 纵向蠕滑率/纵向蠕滑力传递函数系数

表2 横向蠕滑率/横向蠕滑力传递函数系数

表3 自旋蠕滑率/横向蠕滑力传递函数系数

从计算结果中可以看出随着参考纵向、横向或自旋蠕滑率不断增大，对于同一波长蠕滑率的蠕滑力幅值增益不断减小，相位延迟略有增加。传递函数系数随蠕滑率的增加逐渐增加。自旋蠕滑率/横向蠕滑力传递函数常数项和一阶系数随自旋蠕滑率的增加逐渐增加，二阶系数基本不变，三阶系数随自旋蠕滑率的增加逐渐降低。为了对该传递函数进行验证，选取了Lξ=5a的蠕滑率简谐波动工况，蠕滑力对比结果如图3所示。可见使用非稳态传递函数方法和使用Kalker理论计算所得蠕滑力幅值和相位均具有较好的一致性。

图3 Kalker理论与非稳态传递函数方法蠕滑力对比

2.2 其他方向蠕滑率对传递函数的影响

在2.1节中对传递函数进行计算时默认其他方向的蠕滑率为0。线性蠕滑理论公式假设横向蠕滑率对纵向蠕滑力以及纵向蠕滑率对横向蠕滑力均没有影响，但事实上其他方向的蠕滑率对传递函数的性质也会造成一定影响。为研究其他方向蠕滑率对传递函数的影响，设计2个计算工况。

工况1：在纵向参考蠕滑率不变的情况下，施加不同的横向参考蠕滑率(横向波动蠕滑率分量为0)，分别为0.1‰,1.0‰,2.0‰,3.0‰，从而研究横向蠕滑率对纵向蠕滑率/纵向蠕滑力传递函数的影响。

工况2：在横向参考蠕滑率不变的情况下，施加不同的纵向参考蠕滑率(纵向波动蠕滑率分量为0)，从而研究纵向蠕滑率对横向蠕滑率/横向蠕滑力传递函数的影响。

工况1中，纵向蠕滑率简谐波动时纵向蠕滑力随波长比Lξ/a及横向参考蠕滑率变化如图4(a)所示。工况2中，横向蠕滑率简谐波动时横向蠕滑力随波长比Lξ/a及纵向参考蠕滑率变化如图4(b)所示。

其他方向的参考蠕滑率/‰ 0.1 1.0 2.0 3.0图4 蠕滑力随波长比及其他方向参考蠕滑率 变化的Nyquist图

可采用一阶传递函数1/[a0+a1(jw)]对图4进行拟合。工况1采用纵向参考蠕滑率计算出的传递函数系数见表4。工况2采用横向参考蠕滑率计算出的传递函数系数见表5。

表4 工况1传递函数系数

表5 工况2传递函数系数

工况1中，采用多种纵向参考蠕滑率计算的传递函数系数见图5；工况2中，采用多种横向参考蠕滑率计算的传递函数系数见图6。

图5 工况1传递函数系数

图6 工况2传递函数系数

结果表明随着横向参考蠕滑率的增加，同一波长纵向蠕滑率的蠕滑力幅值增益不断减小，相位延迟略有增加。且横向蠕滑率对幅值增益的影响随纵向蠕滑率的增加逐渐减小。当纵向蠕滑率达到3.0‰时，横向蠕滑率在0.1‰,1.0‰,2.0‰,3.0‰时的幅值增益均十分接近。传递函数系数a0和a1均随纵向蠕滑率以及横向蠕滑率的增加而增加。不同纵向参考蠕滑率时横向蠕滑率/横向蠕滑力传递函数规律与其一致。

3 结论

1)相对于传统的稳态滚动接触模型，非稳态滚动接触模型能够体现接触参数高频变化对接触力的影响，更适用于研究高速铁路在高频激励作用下的轮轨瞬态响应问题。使用非稳态滚动接触模型计算蠕滑率简谐激励时蠕滑力幅值增益减小，相位滞后。且随着蠕滑率简谐波动波长比Lξ/a的减小，蠕滑力幅值增益减小程度越大，相位滞后越多。

2)使用非稳态传递函数方法和Kalker三维滚动接触精确解方法计算时，蠕滑力幅值和相位均具有较好的一致性。随着参考蠕滑率不断增大，对于同一波长蠕滑率的蠕滑力幅值增益不断减小，相位延迟略有增加。

3)在纵向蠕滑率简谐激励时，横向参考蠕滑率的增大会减小纵向蠕滑力的幅值增益，但对纵向蠕滑力的相位滞后影响不大。这种幅值增益的减小效果随着横向参考蠕滑率的增大而逐渐减小。在横向蠕滑率简谐激励时，纵向参考蠕滑率对横向蠕滑力的影响与之相似。