脂介质条件下轮轨摩擦因数数值模拟分析

时间：2024-07-28

刘思磊许玉德刘一鸣

1.同济大学道路与交通工程教育部重点实验室，上海201804；2.深圳市综合交通设计研究院有限公司，广东深圳518003

近年来，重载铁路向着大轴重、长编组、高密度方向发展，轮轨状态恶化加剧。运用轮轨摩擦控制剂来控制轮轨间摩擦是延长轮轨使用寿命、降低轮轨间磨耗与损伤的有效手段之一。

早在2004年，就有日本、加拿大、西班牙等国家开始对线路进行润滑试验，并获得了较大的经济效益［1-3］。对于轮轨摩擦因数的变化规律，日本的伴巧、松本耕輔、深貝晋也等［4-6］通过在二圆筒试验机上的试验，发现摩擦因数的改变量与转动周数呈指数关系。近年来，国内学者通过试验进行了相关研究。刘腾飞［7］利用滚动摩擦磨损试验机研究了不同介质对轮轨黏着的影响。李伟［8］发现，进行曲线摩擦控制可改善轨距面的磨耗，增加车辆的曲线通过能力，大大降低轮轨横向力。张信等［9］设计了试验，分析丙烯酸树脂乳液的加入量对涂覆层的影响。宋靖东等［10］分析了固体润滑剂在轮轨间的摩擦性能和有效工作时间情况。王珊珊［11］在点接触试验机上通过球盘接触试验，研究了滑滚比对润滑脂纤维团在接触区分布的影响。

既有文献中关于摩擦控制条件下摩擦因数的影响机理研究较少，且缺少对摩擦因数改善效果的定量研究。本文在调研国内外铁路摩擦控制及摩擦学相关研究的基础上，进行脂介质条件下的轮轨摩擦因数分析与计算。

1 计算模型

轮轨摩擦控制属于弹性流体动压润滑问题，即在流体动压润滑的基础上考虑了弹性变形，而弹性变形因素的比重与膜厚比有关。在求解时，一般将Hertz接触压力分布作为初始压力分布。本文引用弹流理论基本方程，引入脂介质的流变本构模型并考虑接触中热效应等问题。

1.1 脂介质下弹性流体动压润滑理论基本假设

对流体的摩擦状态进行研究时，所用的基本方程是Reynolds方程和膜厚方程［12］。其中Reynolds方程是从黏性流体力学的基本方程出发，作了一定的假设而导出的微分方程。这些基本假设是：①液膜厚度很小，仅有数千纳米，因此沿液膜厚度方向，液膜压力为常数；②轮轨摩擦表面的曲率半径远大于液膜厚度，因此接触面可展为平面；③摩擦控制剂视为脂类，剪应力与剪应变率呈非线性关系；④液膜中流体作层流运动；⑤流体的惯性力忽略不计。

1.2 基本方程

1.2.1 Reynolds方程

Reynolds方程是基于动量方程和连续性方程推导出来的流体动力润滑理论，于1886年由Reynolds提出，是弹流问题中最基本的方程［11］。加入Ostwald模型后，可推导出脂介质的一维Reynolds方程，即

其中，Ostwald模型为

式中：n为流体指数；ρ为流体密度；h为液膜厚度；Φ为塑性黏度；p为流体压力；us为上下接触表面的平均速度；τ为剪应力；γ˙为剪应变率。

入口区的边界条件为p(x0)=0；出口区的边界条件为其中，x0表示入口坐标，xe表示出口坐标。

1.2.2 膜厚方程

膜厚是Reynolds方程的一个重要子项［12］。两个半径分别为R1和R2的圆柱体作线接触时，若两圆柱未变形时的间距为h0，考虑到表面的弹性变形，实际的线接触膜厚公式为

式中：R可由得到；v(x)为弹性变形。

如图1所示，v(x)由弹性力学中表面上的法向分布荷载在半无限大平面上引起的表面变形公式得到，即

式中：p(s)为荷载分布函数；s是x轴上的附加坐标，表示任意线荷载p(s)ds与坐标原点的距离；实际计算中，c为常数，可以同式（3）中的h0合并。

1.2.3 荷载平衡方程

求解弹流问题是在给定荷载条件下进行的，所以求出的压力必须满足荷载平衡条件。当轮轨接触处于部分膜弹流状态时，荷载由固体粗糙峰荷载We和液膜荷载Wa共同承担，荷载方程［13］为

式中：pA为固体A上粗糙峰的压力。

根据Greenwood-Tripp接触理论，pA的计算方程为

当h∕σ在0.000 1～0.010 0之间变化时，kA的取值为0.000 3～0.003 0。由于F5∕2(h∕σ)中积分的计算较复杂，因此采用拟合公式进行计算，即

1.2.4 摩擦因数计算方程

本文弹流问题的荷载由摩擦控制剂的动压力We和液膜荷载Wa共同承担，当求出We和微凸体峰接触所承担的荷载Wa后，总的摩擦因数f的表达式为［14］

式中：μe为脂类摩擦控制剂全膜弹流状态下的摩擦因数，轨顶使用的脂类摩擦控制剂μe约为0.30；μa为微凸体峰接触摩擦因数，一般由试验得到，本文取μa=0.10。

2 数值计算

为方便数值求解，改善方程的适用性，须采用多重网格法求解弹性变形，然后用压力迭代修正。

2.1 多重网格法求解弹性变形

弹流问题的关键是将弹性变形方程叠加到膜厚方程上，弹性变形方程的求解要解偏微分方程，求解过程中选择合适的网格常常是困难的，使用稀疏网格得到的解精度太低，而且对于非线性问题常得不到收敛的解。20世纪60年代初提出了一种求解代数方程组的高效率迭代方法——多重网格法［15］，可以有效地解决上述困难。

利用多重网格法计算弹性变形主要包含以下步骤：

1）压力下传。考虑压力变化，将细网格上的压力插值后传递到粗网格上。

2）积分系数下传。采用相应计算公式将细网格的积分系数下传给粗网格。

3）粗网格积分。

4）积分值插值回代。部分未计算积分值的细网格节点上，如果与粗网格节点重合则作映射运算，如果不重合则作插值运算。

5）细网格修正。

2.2 压力迭代修正

在初始压力分布下（本文为Hertz接触压力）应用多重网格计算方法得到弹性变形后，可以计算膜厚和黏度值，并带入Reynolds方程求解得到新的压力分布，对前一次的压力分布进行迭代修正，然后带入能量方程求解温度。利用新的温度修正黏度值，再迭代求解压力。重复此过程，直至达到精度要求，迭代结束。最终得到所求的压力分布和含弹性变形的膜厚分布。

压力迭代修正同时使用Gauss-Seidel迭代和Jacobi双极子迭代方法［13］。在低压力区，选择Gauss-Seidel迭代法；在高压力区，为避免发散，选择Jacobi双极子迭代法。

1）Gauss-Seidel迭代法

修正时须在本点加上修正量。表达式为

2）Jacobi双极子迭代法

修正时须在本点加上修正量和在前点减去修正量，表达式为

式中：c2为松弛因子。

3 计算结果分析

按照上述计算流程求解基本方程，依据计算结果讨论脂类摩擦控制剂的黏度、列车轴重及列车速度对摩擦因数的影响。

3.1 脂介质黏度对摩擦因数的影响

基于经典弹流理论和脂介质本构方程，计算黏度为10～500倍水黏度时的摩擦因数，模拟中等黏度的脂类摩擦控制剂。

表面粗糙度为6 μm时，选取脂介质黏度分别为0.01、0.10、0.30、0.50 Pa·s（依次对应10、100、300、500倍水黏度）计算得到轮轨接触区压力分布情况及膜厚变化情况见图2。其中，坐标进行了无量纲化处理；横坐标-1.0～1.0为接触区。

由图2可知：

图2 不同脂介质黏度下接触压力和液膜厚度变化曲线

1）从接触区宽度看，随着脂介质黏度增加，液膜和粗糙峰的压力接触区入口逐渐向前移动，导致接触区宽度增加。当脂介质黏度为10倍水黏度（0.01 Pa·s）时，液膜压力和粗糙峰压力的接触区宽度与Hertz接触区宽度基本一致；当脂介质黏度增加至500倍水黏度（0.50 Pa·s）时，液膜压力的接触区宽度约为1.5倍Hertz接触区宽度，粗糙峰压力的接触区宽度约为1.2倍Hertz接触区宽度。

2）从液膜接触压力看，随着黏度增加，在出口区附近，液膜接触压力会出现逐渐明显的二次尖峰。

3）从液膜厚度看，在膜厚出口处均有颈缩现象，随着黏度增加液膜厚度逐渐增加，使得液膜承担的荷载比重上升，粗糙峰承担的荷载比重下降。

不同轮轨表面粗糙度情况下，黏度变化对摩擦因数的影响见图3。

图3 不同表面粗糙度下黏度变化对摩擦因数的影响

由图3可知，随着脂介质黏度增加，摩擦因数呈增加趋势，且在轮轨表面比较光滑时，摩擦因数更早地逼近临界值0.3。产生这一现象的原因是，当黏度增加时，脂介质所形成的液膜将承担更多的荷载，根据式（7），接触区的摩擦因数由液膜承担的荷载比重决定，而脂介质的摩擦因数大于粗糙峰的摩擦因数，因此引起摩擦因数增加。而摩擦因数的上限是全膜弹流状态下的0.3，因此摩擦因数随黏度的增加趋势将趋于平缓。

3.2 轴重对摩擦因数的影响

表面粗糙度为6 μm时，选取轴重分别为10、14、18、22 t，计算得到轮轨接触区压力分布情况及膜厚变化情况见图4。其中，坐标进行了无量纲化处理；横坐标-1.0～1.0为接触区。

图4 不同轴重下接触压力和液膜厚度变化曲线

由图4可知：

1）从接触区宽度看，随着轴重增加，液膜和粗糙峰的压力接触区入口逐渐向后移动，导致接触区宽度缩小。当轴重为10 t时，液膜压力的接触区宽度约为1.9倍Hertz接触区宽度，粗糙峰压力的接触区宽度约为1.4倍Hertz接触区宽度；当轴重为22 t时，液膜压力的接触区宽度与Hertz接触区宽度基本一致，粗糙峰压力的接触区宽度约为1.25倍Hertz接触区宽度。

2）从液膜接触压力看，随着轴重增加，液膜接触压力的二次尖峰逐渐向出口区移动，当轴重大于22 t以后，尖峰不再明显，总应力分布已经接近Hertz应力分布，说明此时摩擦效果接近失效。

3）从液膜厚度看，在膜厚出口处均有颈缩现象，随着轴重增加，液膜厚度逐渐减小，使得液膜承担的荷载比重下降，粗糙峰承担的荷载比重上升。

不同轮轨表面粗糙度情况下，轴重变化对摩擦因数的影响见图5。可知，随着列车轴重增加，摩擦因数呈下降趋势。这是因为随着轴重增加，轮轨间脂介质形成的液膜厚度将被压缩，使得粗糙峰承担的荷载比重增加，而粗糙峰的摩擦因数小于脂介质的摩擦因数，因此引起摩擦因数下降。

图5 不同表面粗糙度下轴重变化对摩擦因数的影响

3.3 速度对摩擦因数的影响

表面粗糙度为6 μm时，选取列车速度分别为60、70、80、90 km∕h，计算得到轮轨接触区压力分布情况及膜厚变化情况见图6。其中，坐标进行了无量纲化处理；横坐标-1.0～1.0为接触区。

图6 不同速度下接触压力和液膜厚度变化曲线

由图6可知：

1）从接触区宽度看，列车速度由60 km∕h增加至90 km∕h，液膜压力的接触区宽度由1倍Hertz接触区宽度增加至约1.08倍Hertz接触区宽度，粗糙峰压力的接触区宽度基本保持在1.25倍Hertz接触区宽度，不随速度变化而改变。

2）从液膜接触压力看，随着列车速度增加，液膜接触压力的二次尖峰位置和峰值大小变化不大。

3）从液膜厚度看，在膜厚出口处均有颈缩现象，随着速度增加，液膜厚度逐渐增加，使得液膜承担的荷载比重上升，粗糙峰承担的荷载比重下降。

不同轮轨表面粗糙度情况下，列车速度变化对摩擦因数的影响见图7。可知，随着列车速度增加，摩擦因数呈增加趋势，但增加的幅度很小，当轮轨表面粗糙度处于一般水平10 μm时，速度由50 km∕h增加至100 km∕h，摩擦因数仅由0.260增加至0.263。

图7 不同表面粗糙度下速度变化对摩擦因数的影响

4 试验验证

为验证模型的可靠性以及数值计算结果的正确性，通过脂介质条件下的滚动接触试验，验证轴重因素对摩擦因数的影响，总结不同轴重下摩擦因数的变化规律。为了能够稳定地将脂介质涂覆在接触面上，试件的转速设定为10 r∕min。按照上述模拟准则，试验所采用的试验力及所对应的模拟轴重见表1，其他主要参数见表2。试验现场如图8所示。

表1 模拟轴重与试验力

表2 试验其他参数

图8 试验现场

对轴重分别为10.0、12.5、15.0、17.5、20.0、22.5 t时试验得出的摩擦因数进行分析。其中轴重15 t时的摩擦因数变化曲线见图9。可知：①涂覆脂介质后，轮轨摩擦因数迅速爬升，之后在0.2附近波动。涂覆阶段摩擦因数有时会出现较大波动的原因是，滚动接触的过程中须不断手工涂刷脂介质，涂刷时涂量不稳定会引起摩擦因数波动。②当停止涂覆后，经过约100 s，摩擦因数降低至干摩擦下的0.11附近。

图9 轴重15 t时摩擦因数变化曲线

将试验参数带入到模型中验证。由于脂类材料未标明该参数，须通过试验测定其取值。在滚动接触试验中，最小试验力为1.2 kN，所产生的接触应力等效于现场的8.6 t轴重，而所用试样的粗糙度仅有0.6 μm。在1.2 kN试验力、10 r∕min的工况下进行了5组试验，5组试验的摩擦因数最大值的均值为0.22，故取脂介质全膜弹流状态时的摩擦因数为0.22。取涂覆稳定阶段的摩擦因数均值作为该轴重下摩擦因数的试验值。得到试验结果和数值计算结果对比，见表3。

表3 试验结果与数值计算结果对比

由表3可知，试验结果的最大相对误差在4%以内，平均相对误差为2.30%，说明理论计算结果与试验结果较吻合。