铁路机车车辆曲线通过是诸多轮轨伤损问题的根源，比如亟待解决的钢轨波磨、钢轨侧磨、轮缘磨耗以及连续型轮轨滚动接触疲劳等[1-2]，其根本原因是曲线通过时轮轨间存在较大的纵向、横向和自旋蠕滑率，轮轨相互作用较直线段更加复杂、激烈。因此，深入理解曲线通过过程中的轮轨动态相互作用是理解和解决上述问题的关键。本文建立了轮轨瞬态滚动接触有限元模型，研究轮轨曲线通过中的高频/瞬态滚动接触行为及短波钢轨波磨的影响。
相比传统的多体动力学模型，有限元模型可求解轮轨间详细的接触应力、应变及接触斑内黏滑分布等，更适于轮轨间接触分析。利用ANSYS商业软件，Telliskivi和Olofsson[3]建立了求解静态磨耗轮轨间单点、两点接触的有限元模型。采用类似静态有限元模型，Sladkowski和Sitarz[4]分析了不同车轮型面和冲角对接触斑以及接触应力的影响。侯传伦[5]建立了三维弹塑性半轮对有限元模型，分析了实际磨耗状态下轮轨的静态接触。需说明，以上有限元模型均仅针对静态轮轨接触，虽可考虑不同横移，但未将曲线钢轨考虑在内。
基于隐-显式结合建模技术的三维轮轨瞬态滚动接触有限元建模方法，在过去20年左右的时间内逐渐发展成熟[6-12]。赵鑫等[7,9,11]开发的高速轮轨三维瞬态半轮对滚动接触模型可充分考虑轮轨表面三维不平顺、不同牵引/制动载荷和低黏着等复杂环境，于时域内模拟直线轨道上速度高至500 km/h下轮轨间瞬态滚动接触行为。近来，刘超等[13]进一步扩展上述模型，建立了可以考虑横向蠕滑和大自旋的全轮对直线段三维瞬态滚动接触模型，并系统分析了单侧钢轨波磨对两侧轮轨瞬态响应的影响。
本文进一步扩展上述全轮对瞬态滚动接触有限元模型至曲线段轨道，将曲线超高、弯曲钢轨、轮对横移及侧滚等因素考虑在内。以R300 m圆曲线轨道为例，详细分析了自由轮对滚过无不平顺的光滑钢轨和单侧钢轨波磨段时两侧轮轨间的瞬态滚滑行为，为深入研究车辆曲线通过行为和曲线段轮轨损伤提供了有效的分析手段。
1 曲线段三维瞬态滚动接触模型
1.1 模型概述
使用ANSYS/LS-DYNA建立了图1所示的考虑R300 m圆曲线轨道和全轮对的三维轮轨瞬态滚动接触有限元模型，涉及参数均取自我国某地铁系统(无砟轨道，最高运行速度120 km/h)。一条自由轮对置于曲线轨道上，可设置不同的轮对横移和侧滚角，以充分模拟轮轨间纵向、横向和自旋蠕滑率及相关的滚滑行为，两簧上质量通过一系悬挂与车轴两端在轴箱位置上相连，曲线超高、弯曲钢轨等均考虑在内。如此，忽略了轮对间相互作用、簧上质量本为一体及曲线通过时的簧载质量转移现象(簧上质量对高频轮轨作用的影响可以忽略[14])。
模型轨道中心线弧长12.6 m或包含20组扣件(扣件间距在中心线上为0.625 m)，超高100 mm，轨底坡设为1∶40。左、右侧一系悬挂各由8组(1×8)平行分布的弹簧和阻尼单元模拟，每一组钢轨扣件由15组(3×5)分布于橡胶垫范围内的弹簧和阻尼单元阵列模拟，每一组弹簧和阻尼单元由空间上重合的一个弹簧和一个阻尼单元构成。真实几何的轮轨采用8节点六面体单元划分，计算时采用减缩积分格式，最小单元尺寸1 mm[7] (图1(b)中细化分析区BC段，采用基于罚函数法的“面-面”接触算法进行轮轨接触计算)。车轮踏面廓形LM，钢轨CN60，采用相同网格离散左、右两侧轮轨，以降低计算误差。该模型总节点数和总单元数分别约为260万和235万，主要参数见表1。
为表述方便，建立图1所示的笛卡尔坐标系oxyz，坐标原点位于轮对初始位置A处，即轨道中心在过两轨底中心线平面上的投影，x轴沿初始滚动方向。轨道假设为右曲线，面向运行方向的左侧为外轮/轨，右侧为内轮/轨。
图1 曲线段上的三维轮轨瞬态滚动接触模型 /m
Fig.1 An overview and a zoom-in of mesh of 3D transient wheelset-rail rolling contract model on curve
表1 模型参数表
Table 1 Parameters of the FE model
参数 数值簧上总质量m/kg 12857.2轮对质量m/kg 1142.8车辆一系悬挂(单侧) 刚度系数K/(MN/m) 0.98阻尼系数C/(kN·s/m) 8.2轮对、钢轨 杨氏模量E/GPa 205.9密度ρ/(kg/m3) 7790泊松比 0.3阻尼常数 1.0×10-4橡胶垫 刚度系数K/(MN/m) 40阻尼系数C/(kN·s/m) 50轨道板 杨氏模量E/GPa 32.5密度ρ/(kg/m3) 2400泊松比 0.24摩擦系数 0.5
1.2 典型模拟过程
模拟过程与文献[9]类似，首先采用隐式时间积分求得车轮在初始位置A(图1(b))的静态位移场，然后，以此位移场初始化轮轨应力状态，设置其他初始条件和载荷，再通过显式时间积分求解轮对以指定速度和牵引系数从A沿曲线轨道滚向C的瞬态滚动接触行为。与文献[9]中模型相比，除了曲线轨道之外，本模型还有如下差异或改进：1) 轮对初始横移和侧滚角不再为零，而是施加了自由轮对稳态通过曲线时的结果，轮对初始冲角依然为零(模拟非导向轮对)；2) 初始静态位移场计算时，也施加了相应的横移及侧滚角。瞬态计算的边界条件如下：一系悬挂的上、下节点(对应簧上质量与车轴轴心节点)在横、纵向耦合，保证簧上质量与轮对一起运动；钢轨和轨道板纵向两端对称约束，以模拟其无限长特性；轨道板两侧面横向约束，底面全约束；一系悬挂和钢轨扣件只保留垂向自由度。
1.3 轮对横移和侧滚角
本文仅针对自由轮对稳态通过圆曲线段工况，轮对横移量[15] (y0，轮对中心距线路中心线的距离)的公式如下。由轮对在纯滚线上的接触几何(见图2)可得到式(1)，联立式(2)中踏面等效锥度的定义，可得轮对横移计算式(3)。
式中：r0为车轮名义滚动半径；rl(r)为左(右)轮的实际滚动圆半径；R为曲线半径；2b为轮对左右名义滚动圆横向跨距；λ为等效锥度。
将式(2)和式(3)分别绘制于图3，两线的交点即为轮对稳态通过圆曲线段时的横移量，可得，本文LM车轮与CN60钢轨匹配下轮对稳态通过R300 m圆曲线时的横移量为5.3 mm。求得轮对横移后，其相应侧滚角采用由“试凑逼近法”编写的轮轨刚性接触点对程序求得。图4显示了LM车轮与CN60钢轨匹配时侧滚角与横移量的关系，可知横移5.3 mm对应的侧滚角为0.0359º。
图2 自由轮对纯滚过曲线的几何示意图
Fig.2 Geometric diagram of wheelset passing over curves
图3 LM车轮与CN60钢轨匹配下轮对稳态通过R300 m圆曲线段时横移量的确定方法
Fig.3 Method for determining the lateral shift of the LM wheel marched with CN60 rail when the wheelset passes over the R300 m circular curve
图4 LM车轮与CN60钢轨匹配时侧滚角与横移量的对应关系
Fig.4 Correspondence between roll angle and lateral shift when LM wheel is matched with CN60 rail
1.4 波磨模拟
网格划分完成后的轮轨表面是无几何不平顺的光滑轮轨，三维波磨通过修改BC段表面相关节点的坐标施加。鉴于所模拟地铁线路上某R300 m曲线段的内、外轨均存在40 mm短波波磨，如图5所示。在这一实测结果的基础上，将该波磨抽象和一般化为波长40 mm的余弦形波磨，波深取0.03 mm，且只施加于一股钢轨(内轨或外轨)，具体施加的波磨几何见图6。需说明，波磨深度沿横向呈抛物线分布，波宽为60 mm。
2 光滑轮轨结果
2.1 稳态曲线通过
轮对以50 km/h滚过R300 m光滑圆曲线段(均衡速度为50.42 km/h)，轮对横移量和侧滚角取上述5.3 mm和0.0359º，牵引系数[16]取0.3，轮轨力结果见图7。可以看出，两侧垂向轮轨力在BC段均基本稳定在静轮重70 kN上下；纵向轮轨力稳定在牵引系数对应的20 kN上下，且外侧轮轨大于内侧轮轨，以提供曲线通过所需的扭矩；内轨横向力很小，稳定在零附近，外轨上由于接触角及其导致的自旋等的存在，其值约为4.5 kN，也基本稳定。
图5 我国某地铁线路曲线段(R300 m)上内、外轨波磨的CAT实测结果
Fig.5 Corrugation on the inner and outer rails measured by CAT on a curve of a subway in China (R300 m)
图6 模型中施加的短波波磨几何
Fig.6 Geometry of short-pitch corrugation applied in model
图7 50 km/h通过R300 m曲线时的轮轨力结果(光滑轮轨)
Fig.7 Wheel and rail contact forces when passing by a R300 m curve at 50 km/h (smooth rail)
图8对比了t=0.0527 s时刻内、外轨接触斑的黏滑分布。可以发现，外轨滑移区明显大于内轨，且更靠近轮缘(即轨距角，外轨轨距测量点y=705.4 mm)，这是因为外轨侧更大的纵向力和横向力(见图7)。相应的接触斑内微滑分布示于图9，其中箭头方向和长度分别表征相对滑移速度的方向和大小，从中可以看出由非零接触角导致的自旋蠕滑，且微滑分布与图8中粘滑分布相吻合。
图10给出了上述工况中施加于内、外侧钢轨轨头的摩擦功分布。若假设磨耗与摩擦功成正比，该图也显示了两轨磨耗分布情况。发现，外轨磨耗远大于内轨，最大值约为内轨的3.1倍，且磨耗位置比内轨更靠近轨距角，与现场观测一致。
图8 50 km/h通过R300 m曲线过程中t=0.0527 s时刻的接触斑内黏滑分布(光滑轮轨)
Fig.8 Stick-slip distributions in the contact patch at t=0.0527 s when the wheelset rolls over a R300 m curve at 50 km/h(smooth rails)
图9 50 km/h通过R300 m曲线过程中t=0.0527 s时刻的接触斑内微滑分布(光滑轮轨)
Fig.9 The micro-slip distributions in the contact patch at t=0.0527 s when the wheelset rolls over a R300 m curve at 50 km/h (smooth rails)
图10 50 km/h通过R300 m曲线过程中轨头的摩擦功分布(光滑轮轨)
Fig.10 The frictional work on the rail head when the wheelset rolls over a R300 m curve at 50 km/h (smooth rails)
需指出，上述结果均符合稳态曲线通过的基本特征，说明所建立的模型能合理地模拟轮轨曲线通过行为。
2.2 横移影响
仍以上述R300 m圆曲线为例，由线路激扰导致的不同横移和相应侧滚角时轮轨力在结果分析段内的平均值如图11所示。除初始横移量及相应侧滚角之外，其他参数均保持不变，详见表2。由图11可见，对于垂向轮轨力，随着轮对横移量的增大，外轨侧增载，内轨侧减载，在稳态横移5.3 mm时两者相同；纵向轮轨力在横移0 mm时，内、外轨差异最小，且随着横移增大而单调递增，即轮对转向力矩逐渐增大。内、外侧轮轨横向接触力差异随轮对横移量的增大也愈加明显，且外轨横向接触力幅值明显增加，这是由于外侧轮轨接触点对逐渐向轮缘根部移动，横向和自旋蠕滑率相应增大所致。
图12展示了t=0.0527 s时刻4种横移下轮轨接触斑形状和黏滑分布结果。不难看出，随着横移的增加，外轨接触斑明显向轨道内侧(轨距角)移动，且接触斑面积先增大后减小；内轨接触斑向轨道外侧(轨顶)扩展，接触斑面积单调递增。
表2 不同横移工况计算参数
Table 2 Calculation parameters for different traverse condition
工况序号 轮对横移/mm 轮对侧滚/(º) 速度/(km/h)1 0 0 50 2 3 0.0181 50 3 5.3 0.0359 50 4 8 0.0705 50
图11 轮轨力随轮对横移的变化
Fig.11 Variation of forces between wheel and rail with different lateral shifts
图12 不同横移工况中t=0.0527 s时刻的接触斑形状、面积(图中用数字标出)及黏滑分布
Fig.12 The shape, areas and stick-slip distributions at t=0.0527 s with different lateral shifts
由图13所示的不同横移下同一时刻接触应力沿纵向路径Y (定义为经过各自接触斑内最大接触应力点的钢轨表面纵向路径)的2D分布可见，外轨法、切向接触应力均先减小后增大，与接触斑面积变化趋势相对应，且横移8 mm时，所选取纵向路径上的切向接触应力接近饱和。内轨法、切向接触应力随横移增加而逐渐减小，这是总法向、切向力减小而接触斑面积不断增大的结果。需说明，为方便对比，图13只展示了二维结果，并将接触斑起点平移至同一纵向位置(不同横移时相同时刻下接触斑的位置稍有差异)，且法向应力结果为法向接触应力与摩擦系数0.5的乘积(即极限切向接触应力)。
图14给出了细化分析区内轨头纵向路径L (经过轨头摩擦功最大值)上的摩擦功平均值随横移的变化。该结果显示，外轨摩擦功在横移3 mm时出现最小值9.87 mJ，随着横移的进一步增加，摩擦功显著增加，例如横移8 mm时的摩擦功高达56.94 mJ(对应着轨距角处的严重磨耗)，而内轨摩擦功随横移的增大而单调递减。需说明，不同横移工况所对应的纵向路径L的横向位置不同。
图13 不同横移下t=0.0527 s时刻轮轨接触应力沿纵向路径Y的分布
Fig.13 The contact stresses between wheel & rail along longitudinal path Y at t=0.0527 s with different lateral shifts
图14 内轨、外轨的轨头摩擦功随轮对横移的变化
Fig.14 Variation of the frictional work on the rail head with different lateral shifts
综上，当轮对横移因激扰超过稳态横移时，外轨的接触力、接触应力、摩擦功均急剧增大。以此，为避免外轨发生严重磨耗及损伤，应尽量减少车辆曲线通过时的横向激扰，以降低大横移的出现频率。
3 短波波磨结果
文献[13]的结果显示，直线轨道上一侧钢轨波磨激励的高频轮轨相互作用可传递至另一侧，引发另一侧钢轨的不均匀磨损。在本文模型上于内轨或外轨分别施加图6所示的波长40 mm的钢轨波磨，计算分析了不同运行速度下曲线段单侧短波波磨对两侧轨头磨损的影响。
当短波波磨存在于外轨时，轮对以70 km/h速度通过时内、外轨的轨头摩擦功分布如图15所示。可以发现，外轨明显的不均匀摩擦功会导致内轨(无波磨侧)出现幅值较小的不均匀摩擦功(约为外轨的1/12)。若假设不均匀摩擦功必然导致不均均磨损，外轨短波波磨会引发内轨逐渐出现类似波长的波磨。
短波波磨分别施加于外轨或内轨上时，轮对以不同速度通过上述R300 m圆曲线所造成的轨头摩擦功波动幅值如图16和图17所示。此波动仅考虑与波磨通过频率相关的的波动成分。由图16可见，当短波波磨仅发生在外轨时，外轨摩擦功波动幅值的极小值发生在40 km/h，并非在近似均衡速度50 km/h，70 km/h速度下不均匀摩擦功的波动幅值达到最大值8.05 mJ，是40 km/h时的2.6倍。相对而言，内轨不均匀摩擦功波动幅值在速度40 km/h~60 km/h范围内变化不大，在30 km/h和70 km/h时的数值分别为0.63 mJ和0.69 mJ，但依然仅为外轨最大值的8%左右，或可导致内轨出现类似波长的轻微钢轨波磨。
图17为短波波磨存在于内轨时的相应结果。不难发现，内、外轨摩擦功波动幅值规律相似，且均在40 km/h时最小，同时无波磨侧的外轨摩擦功最大波动幅值在70 km/h时达到最大，是外轨波磨时内轨最大波动幅值的1.9倍。这说明，短波波磨存在于内侧钢轨时，比存在于外轨，可使所激励的振动更有效地传至无波磨侧钢轨，最终导致双侧波磨现象。另外，内轨波磨造成的内轨侧最大摩擦功波动幅值远低于外轨波磨造成的外轨侧最大值。造成这些现象的原因是外轨侧更高的蠕滑力(见图9)。
综合图16和图17的结果可以发现，速度对两侧不均匀磨损有不可忽略的影响。从两个结果来看，两侧摩擦功波动幅值均在40 km/h~50 km/h间某速度下达到最小值，最有利于降低波磨的发展速率。
图15 轮对以70 km/h通过外轨存在短波波磨的R300 m圆曲线时的轨头摩擦功分布
Fig.15 The frictional work on the rail head during curving at 70 km/h (with short-wave corrugation on outer rail, R300 m)
图16 外轨短波波磨造成的轨头摩擦功波动幅值随速度的变化
Fig.16 Amplitudes of the frictional work fluctuations on rail heads at different speeds when short-pitch corrugation exists on the outer rail
图17 内轨短波波磨造成的轨头摩擦功波动幅值随速度的变化
Fig.17 Amplitudes of the frictional work fluctuations on rail heads at different speeds when short-pitch corrugation exists on the inner rail
4 讨论
众所周知，曲线通过会造成不可忽略的簧载转移现象，继而对轮轨接触产生重要影响。如何将与超高、速度等参数相关的簧载转移合理地考虑在内，是下一步研究中应该考虑的一个重要问题。初步计算表明，外轨增载而内轨相应减载时，内、外轨之间的磨耗差异会较本文所示结果更大，即外侧轮轨的服役环境较文中模拟工况更加恶劣，更多研究正在进行中。另外，考虑到现场中波磨波长随速度的变化等现象，未来也应开展更系统的速度影响研究。
5 结论
建立了考虑曲线段的全轮对三维轮轨瞬态滚动接触有限元模型，详细分析了无几何不平顺钢轨和单侧钢轨波磨两种工况下的自由轮对曲线通过行为。基于所示瞬态滚动接触结果，得出如下结论：
(1) 该模型可于时域内模拟真实几何自由轮对通过曲线时的轮轨瞬态滚动接触行为。
(2) 以50 km/h通过光滑钢轨曲线时(均衡速度50.42 km/h)，外轨磨耗远大于内轨，最大磨耗值约为内轨的3.1倍，且磨耗位置比内轨更靠近轨距角。
(3) 当轮对横移因激扰超过稳态横移量时，外轨的轮轨力、接触应力和摩擦功会急剧增大。
(4) 曲线段单侧钢轨存在波长40 mm短波波磨时，也会引发另一侧无波磨钢轨的不均匀磨损，且存在于内轨时，其引起的另一侧不均匀磨损的幅值更大(差别可达1.9倍)。速度对两侧钢轨不均匀磨损有不可忽略的作用，就本文所模拟工况而言，40 km/h~50 km/h间某曲线通过速度(略低于均衡速度)下波磨发展速率最低。
[1]王开云, 翟婉明, 刘建新, 等.提速列车与曲线轨道的横向相互动力作用研究[J].中国铁道科学, 2006,26(6): 38-43.Wang Kaiyun, Zhai Wanming, Liu Jianxin, et al.Research on the lateral dynamic interaction between speed-increased train and curve track [J].China Railway Science, 2006, 26(6): 38-43.(in Chinese)
[2]熊嘉阳, 金学松.铁路曲线钢轨横向凹坑对初始波磨形成的影响[J].工程力学, 2006, 23(06):135-141.Xiong Jiayang, Jin Xuesong.Effects of lateral dent on curved rail on rail corrugation [J].Engineering Mechanics, 2006, 23(06): 135-141.(in Chinese)
[3]Telliskivi T, Olofsson U.Contact mechanics analysis of measured wheel-rail profiles using the finite element method [J].Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail and Rapid Transit,2001, 215(2): 65-72.
[4]Sladkowski A, Sitarz M.Analysis of wheel-rail interaction using FE software [J].Wear, 2005, 258(7):1217-1223.
[5]侯传伦.重载铁路曲线段磨耗状态下轮轨相互作用分析[D].成都: 西南交通大学, 2009.Hong Chuanlun.Analysis of worn wheel/rail interaction on heavy-haul railway curve [D].Chengdu: Southwest Jiaotong University, 2009.(in Chinese)
[6]常崇义.有限元轮轨滚动接触理论及其应用研究[D].北京: 中国铁道科学研究院, 2010.Chang Chongyi.A study on wheel/rail rolling contact theory based on finite element method and its applying[D].Beijing: China Academy of Railway Sciences, 2010.(in Chinese)
[7]Zhao X.Dynamic wheel/rail rolling contact at singular defects with application to squats [D].Delft: Delft University of Technology, 2012.
[8]肖乾.轮轨滚动接触弹塑性分析及疲劳损伤研究[D].北京: 中国铁道科学研究院, 2012.Xiao Qian.The elasto-plastic analysis and fatigue damage research of wheel/rail rolling contact [D].Beijing: China academy of railway sciences, 2012.(in Chinese)
[9]赵鑫, 温泽峰, 王衡禹, 等.三维高速轮轨瞬态滚动接触有限元模型及其应用[J].机械工程学报, 2013,49(18): 1-7.Zhao Xin, Wen Zefeng, Wang Hengyu, et al.3D transient finite element model for high-speed wheel-rail rolling contact and its application [J].Chinese Journal of Mechanical Engineering, 2013, 49(18): 1-7.(in Chinese)
[10]宋华, 杨建, 张月, 等.非线性稳态曲线通过时轮轨滚动接触的数值求解方法[J].中国铁道科学, 2015,36(5): 80-85.Song Hua, Yang Jian, Zhang Yue, et al.Numerical solution method for wheel-rail rolling contact of nonlinear steady-state curve negotiating [J].China Railway Science, 2015, 36(5): 80-88.(in Chinese)
[11]Zhao X, Zhao X, Liu C, et al.A study on dynamic stress intensity factors of rail cracks at high speeds by a 3D explicit finite element model of rolling contact [J].Wear,2016, 366-367(11): 60-70.
[12]赵小罡, 赵鑫, 温泽峰, 等.轮轨黏着系数对钢轨直裂纹瞬态扩展行为的影响[J].工程力学, 2018, 35(5):248-254.Zhao Xiaogang, Zhao Xin, Wen Zefeng, et al.Influence of wheel-rail adhesion coefficient on transient propagation of a vertical rail crack [J].Engineering Mechanics, 2018, 35(5): 248-254.(in Chinese)
[13]刘超, 赵鑫, 赵小罡, 等.单侧钢轨波磨对两侧轮轨瞬态响应的影响分析[J].机械工程学报, 2017, 53(22):117-124.Liu Chao, Zhao Xin, Zhao Xiaogang, et al.Analyses of transient wheel-rail interactions excited by unilateral rail corrugation [J].Chinese Journal of Mechanical Engineering, 2017, 53(22): 117-124.(in Chinese)
[14]Chaar N, Berg M.Simulation of vehicle-track interaction with flexible wheelsets, moving track models and field tests [J].Vehicle System Dynamics, 2006, 44(Supp1):921-931.
[15]罗仁, 石怀龙.铁道车辆系统动力学及应用[M].成都:西南交通大学出版社, 2018.Luo Ren, Shi Huailong.Dynamics of railway vehicle systems and application [M].Chengdu: Southwest Jiaotong University Press, 2018.(in Chinese)
[16]刘超, 赵鑫, 安博洋, 等.钢轨短波长波磨处的高速滚动接触分析[J].润滑与密封, 2015, 40(8): 40-46.Liu Chao, Zhao Xin, An Boyang, et al.Study on high-speed wheel-rail rolling contract on short-pitch rail corrugation [J].Lubrication Engineering, 2015, 40(8):40-46.(in Chinese)