轨道交通行业的持续发展极大地缩短了时空距离，并有力带动了世界各地的经济增长。然而，在轨道交通运营里程不断增加的同时，列车运行诱发的轮轨异常振动和噪声问题也越来越受到关注^[1]。为解决这一问题，不少地铁线路引入了减振扣件，如北京地铁多条线路的多个区段采用剪切型减振器扣件降低地铁轮轨异常振动对周围环境的影响^[2-4]，南京地铁1号线换铺嵌套型减振扣件以降低车站及周边建筑的振动噪声^[5]，成都地铁1号线采用GJ-Ⅲ型减振扣件来控制轨道交通噪声和振动^[6]等。虽然减振扣件在一定程度上缓解了异常振动和噪声问题，但会导致减振扣件所在区段产生较为严重的钢轨波磨问题，尤其是科隆蛋扣件(Cologne egg fastener)^[7-9]。科隆蛋扣件属于高弹性减振扣件，垂向静刚度约12 MN/m，在铺设科隆蛋扣件的北京、上海等地地铁线路上均出现了大范围钢轨波磨现象，包括直线和曲线区间；而使用普通扣件(如DTVI2扣件)的直线轨道，绝大多数不会发生波磨，并且当相关部门使用DTVI2扣件更换科隆蛋扣件后，原来波磨严重的直线区间再也没有出现波磨^[10]。钢轨波磨的存在加剧了轮轨系统的高频振动，易导致车辆和轨道部件过早疲劳失效，从而影响列车安全运行^[1]。

关于钢轨波磨的研究，迄今为止已有100余年。文[11-12]根据波长固定机理和损伤机理，将钢轨波磨分为pinned-pinned共振、车辙、其他P2共振、重载、轻轨和特定轨道型式6种类型，并分别给出了相应的治理对策。文[13-14]对20世纪世界范围内铁路线路上的钢轨波磨现象进行了回顾和总结，并提供了丰富的工程信息和理论见解。金学松等^[15]阐述了重载铁路、高速铁路、地铁和城市轨道的钢轨波磨特征及对应的线路结构特征，并将钢轨波磨归因于自激振动理论和反馈振动理论。自激振动理论^[16-18]认为钢轨波磨主要与车辆/轨道系统的固有特性有关，而反馈振动理论^[19-22]则认为钢轨波磨是由轮轨初始不平顺引起，并呈现“不平顺-振动-不平顺”的循环过程^[1]。Sun等^[23-24]从轮轨黏滑振荡过程研究了钢轨波磨的形成原因，结果发现轮轨黏滑过程的频率由一个基频组成，该基频与轨枕通过频率及轮对的扭转和弯曲组合频率相匹配。文[25-26]也从轮轨接触黏滑角度探究了钢轨波磨的形成机制，包括直线线路和曲线线路上的钢轨波磨类型，并结合现场和室内试验对所得结论进行了验证。

由于钢轨波磨发生在轮轨表面接触区域，与微观层面的材料磨耗密切相关，在微观视角下从轮轨接触黏滑特性层面研究钢轨波磨的形成机理更直观且具有说服力。此外，现有文献关于钢轨波磨黏滑振动的计算大多基于半空间假设和稳态理论，对处理轮轨动态相互作用引起的非稳态过程有很大局限性^[27]，且未考虑结构在高频范围内的柔度特性^[1]。有限元法可以在考虑轮轨系统动力学和接触力学的情况下表征实际的轮轨黏滑运动过程^[1]，从而有益于阐释钢轨波磨的形成机理。鉴于此，本文尝试利用有限元法，从轮轨瞬态接触黏滑振动入手分析钢轨波磨的形成过程，以期为钢轨波磨的形成机理提供新见解。

1 波磨现场调研

现场调研区段位于上海某地铁线路，该线路运营车辆为地铁A型车，实测波磨区间为直线区间，采用科隆蛋扣件固定钢轨并实现减振目的，轨下基础为整体轨道板结构。科隆蛋扣件主要利用承轨板与底座之间的硫化橡胶圈的剪切变形获得弹性。由于科隆蛋扣件节点垂向刚度较低，在荷载作用下钢轨的变形较大，因此该扣件降低了轨道结构的整体刚度和轮轨系统的动态冲击。

波磨现场测试采用CAT(corrugation analysis trolley)波磨采集仪，分别测试波磨区间的左右轨。该设备采用惯性基准原理进行测量，探头紧贴轨面进行数据采集，测量重复性不低于0.001 mm，测量分辨率高于0.001 mm，数据采样间隔为2 mm，设备测量速度为0.8~1.2 m/s。钢轨波磨现场照片如图 1所示。

实测轨面不平顺如图 2a所示，其对应的1/3倍频程波长谱如图 2b所示。

图 2中粗糙度级可定义为

$ L=10 \times \lg \left(\frac{r_{\mathrm{rms}}}{r_0}\right)^2 . $

(1)

其中：粗糙度级L由轨面不平顺对应的$ \frac{1}{3}$倍频程波长谱值来表征，dB；r₀为轨面短波不平顺参考值，取值1 μm；r_rms为图 2a中不平顺波深数据进行均方根(root mean square)处理后所得值，即有效值，μm。

由图 2b可知，实测波磨区间波长40~50 mm对应的钢轨粗糙度级超出限值，因此，钢轨波磨的特征波长主要为40~50 mm。该区间车辆运行速度约为55 km/h，根据频率计算公式，可得实测的波磨通过频率为306~382 Hz，波磨通过频率计算公式可表示为

$ f=\frac{1\;000 v}{3.6 \lambda} . $

(2)

其中：v为车辆运行速度，km/h；λ为波磨特征波长，mm。

2 三维轮轨滚动接触有限元模型 2.1 有限元模型的建立

利用有限元软件ABAQUS，参照实际线路情况，建立了采用科隆蛋扣件轨道的三维轮轨滚动接触有限元模型。该模型能够考虑车轮/轨道结构高频范围内的柔度特性、真实接触几何、黏滑振动和瞬态冲击等因素，并使用接地弹簧模拟地基的支撑作用^{[1, 28]}。由于轮对-轨道系统沿线路中心具有对称性，因此仅取该系统的1/2进行建模，有限元模型主要组成部件和参数如表 1所示^[1]。

为详细分析车轮运行过程中的接触状态，主要针对在有无轨面缺陷条件下的轮轨接触黏滑特性进行研究，并在建立模型时考虑无磨耗车轮在直线区间内光滑或含有凹坑缺陷钢轨上的滚动接触行为。车轮初始位置位于科隆蛋扣件上方，无侧滚和冲角，并与钢轨顶面对中接触^[1]。车轮型面为LM，钢轨型面为CHN60，轨底坡为1/40，轨距为1 435 mm；扣件区域纵向长度为200 mm，横向宽度与钢轨轨底宽度一致；轨道板断面尺寸为1 350 mm×300 mm。

轮轨空间运动和材料变形采用Lagrange方法进行描述，车轮、钢轨和轨道板实体结构均采用C3D8R单元进行离散，且假设结构材料为线弹性^[1]。一系悬挂、科隆蛋扣件和接地弹簧均采用弹簧和阻尼单元进行模拟，其中一系悬挂由1组弹簧-阻尼单元连接簧上质量和车轮耦合中心，科隆蛋扣件3个方向分别由300组弹簧-阻尼单元沿扣件区域均匀分布(纵向20组、横向15组)，接地弹簧均匀分布在轨道板底部以连接轨道板和地基，共计75组(纵向25组、横向3组)^[1]。有限元模型整体和局部图如图 3所示。整个车轮耦合约束于车轮中心，并对其施加平动和转动速度，对簧上质量施加相同的平动速度以模拟轮轨滚动接触运动。

三维轮轨滚动接触有限元模型长度为12 000 mm，经测，该长度可有效避免边界条件对模型求解精度的影响^[29-30]。采用隐式+隐式算法求解轮轨滚动接触问题：首先，设置静态接触区长度为75 mm，采用隐式算法求解轮轨初始静态时的接触斑形状和接触压力^[1](施加重力载荷，加载方式为瞬态加载)；其次，考虑到将静态接触解转换至动态接触再次进行隐式求解会引入较大的初始激扰，因此在静态接触区之后设置了长度为1 762.5 mm的动态松弛区，以使初始激扰引起的轮轨系统振动达到稳定；最后，在动态松弛区之后设置了长度为600 mm的动态求解区。需要说明的是，第一次采用隐式算法获取静态接触解是为了使第二次采用隐式算法求解时轮轨系统达到稳定所需的动态松弛区长度减少，从而减少计算量，提高计算效率。静态接触区、动态松弛区和动态求解区在钢轨长度方向上的分布情况如图 3a所示。

为满足模型规模与计算精度和效率之间的最大平衡，静态接触区和动态求解区内的轮轨接触单元尺寸设定为1 mm×1 mm，其他区域进行过渡稀疏网格划分^[1]，完成划分的科隆蛋扣件轨道三维轮轨滚动接触有限元模型共计117.3×10⁴个单元和145.6×10⁴个节点。在模型中建立Descartes坐标系以表示轮轨滚动过程中的相对位移，如图 3a和3b所示。采用隐式算法求解轮轨滚动接触问题时有限元模型的边界条件具体如下：钢轨和轨道板的纵向两端面沿z向对称约束，车轮端面沿x向约束，车轮中心面沿x向对称约束，簧上质量点保留y向和z向的自由度^[1]。

轮轨接触采用“表面-表面”接触算法，滑移公式使用有限滑移以防止主从接触面相对滑移量过大(对于模型最小单元特征尺寸1 mm而言，相对滑移量不大于0.2 mm)，从而减少计算误差^[1]。摩擦行为采用Coulomb摩擦模型进行描述，切向接触使用罚函数法，摩擦系数设置为0.35^[31]；法向接触采用“硬”接触，并允许主从接触面接触后分离以模拟滚动接触工况。通过设定车轮的平动速度和转动速度以及簧上质量点的平动速度，可实现轮轨滚动接触的仿真。

由于瞬时完成速度的施加会促使轮轨系统产生较大的动力学相互作用，为减小上述影响，采用正弦加载函数曲线完成速度的施加，如图 4所示。正弦加载函数可表示为：

$ \boldsymbol{V}_{\mathrm{T}}(t)= \begin{cases}\frac{\boldsymbol{V}_0}{2}\left(1-\cos \left(\frac{{\rm{ \mathsf{ π}}} t}{t_0}\right)\right), & t<t_0 ; \\ \boldsymbol{V}_0, & t \geqslant t_0 ;\end{cases} $

(3)

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{\mathrm{T}}(t)= \begin{cases}\frac{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_0}{2}\left(1-\cos \left(\frac{{\rm{ \mathsf{ π}}} t}{t_0}\right)\right), & t<t_0 ; \\ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_0, & t \geqslant t_0 .\end{cases} $

(4)

其中：V_T(t)和Φ_T(t)分别表示关于时间t的平动速度和转动速度的函数，单位分别为km/h和rad/s，且满足V_T(t)= 0.003 6Φ_T(t)R，R为车轮名义滚动圆半径，取值420 mm；V₀和Φ₀分别为最大平动速度和最大转动速度，分别取值55 km/h和36.38 rad/s；t₀为平动速度和转动速度达到最大值对应的起始时刻，取值0.013 s。

2.2 有限元模型的验证

设定车轮横移量为0，牵引系数(轮轨切向力与轮轨法向力之比^[28])为0，计算车轮以55 km/h的平动速度滚过光滑钢轨时的轮轨接触力，结果如图 5所示。由图 5分析可知，经过一段时间后轮轨接触力均趋于稳定，其中轮轨垂向力稳定在52 286 N左右，轮轨纵向力和横向力基本稳定在0。结果表明：车轮在运行到没有横移和没有激励的直线区间时，经过1 800.0 mm((1 762.5+75.0)/2=1 800.0)的动态松弛后能够基本进入稳定滚动状态，说明本文建立的三维轮轨滚动接触有限元模型计算结果符合实际运行情况，可应用于下一步数值分析。

基于建立的有限元模型开展轮轨接触黏滑特性、轮轨系统的固有振动特性和钢轨纵向磨耗特征分析，具体计算流程如图 6所示。

3 轮轨接触黏滑特性

在三维轮轨滚动接触有限元模型中，通过设置摩擦力阈值，可以区分接触区域的黏着区和滑移区。在黏着区域，节点摩擦力F_i可表示为

$ \boldsymbol{F}_i=k_0 \Delta \boldsymbol{\gamma}_i. $

(5)

其中：k₀为黏着刚度；γ_i为节点相对滑移位移，Δγ_i为节点相对滑移位移的变化量, i取1和2时，分别表示纵向和横向。在滑移区域，节点摩擦力F_i可表示为：

$ \boldsymbol{F}_i=\boldsymbol{F}_{\mathrm{t}} \frac{\Delta \boldsymbol{\gamma}_i}{\Delta \boldsymbol{\gamma}}. $

(6)

$ \boldsymbol{F}_{\mathrm{t}}=\mu \boldsymbol{p}, $

(7)

$ \boldsymbol{F}_{\mathrm{f}}=\sqrt{\boldsymbol{F}_1^2+\boldsymbol{F}_2^2} . $

(8)

其中：F_t为摩擦力阈值；γ为节点相对滑移位移总和，Δγ为节点相对滑移位移总和的变化量；μ为摩擦系数；p为节点法向接触应力；F_f为节点摩擦力合力。当F_f=0时，接触点处于未接触状态；当0＜F_f＜F_t时，接触点处于黏着状态；当F_f≥F_t时，接触点状态则由黏着转变为滑移。下文将主要分析轨面光滑和轨面缺陷工况下的轮轨接触黏滑特性，以期从微观层面解释钢轨波磨的形成机理。

3.1 轨面光滑工况下轮轨接触黏滑特性

当车轮以55 km/h的平动速度运行在光滑轨面上时，动态求解区的轮轨滚动接触黏滑分布图如图 7所示，其中黏合状态表示从接触面绑定于主接触面。需要说明的是，为探究轮轨接触黏滑特性与钢轨波磨的关系，以钢轨表面的接触斑为例进行展示，并通过均匀划分动态求解区的5个代表性时刻(分别对应距离轮轨接触初始位置1.80、1.95、2.10、2.25和2.40 m时的状态)的黏滑分布图予以说明。

由图 7可知，在动态求解区，当车轮平稳滚过光滑钢轨时，轮轨接触斑区域黏滑分布图大致相似，均表现为接触斑前缘被黏着区域占据，接触斑中后缘被滑移区域占据，这与采用稳态动力学计算所得的结果相近^[32-33]，表明本文所建立的有限元模型有效。此外，相较于稳态动力学计算结果，采用瞬态有限元法计算所得的结果更全面地考虑了真实的轮轨接触几何与轮轨结构的柔度特性，因此，其得到的轮轨滚动接触黏滑分布图更接近实际情况。同时，图 7所示5个时刻的黏滑分布均表明轮轨切向力在宏观上未达到饱和，即轮轨接触未达到全滑移状态，因此，车轮一直处于稳定滚动状态，且轮轨系统整体处于稳定状态(如图 5中1 800~2 400 m处图像所示)，这说明在轨面光滑工况下，该轮轨系统不易发生失稳状态，从而不易形成钢轨波磨。

3.2 轨面缺陷工况下轮轨接触黏滑特性

根据3.1节分析可知，在轨面光滑工况下，钢轨表面难以形成波磨，但是在实际线路中却出现了较为严重的钢轨波磨。根据黏滑振动理论^[34-37]可知，轨面固定缺陷能够激发轮轨系统的固有振动特性并赋予钢轨波磨固定波长属性，因此，为进一步从接触黏滑角度探究钢轨波磨的形成机理，需要分析在轨面缺陷工况下的轮轨接触黏滑特性。

设定钢轨表面含有凹坑缺陷，位于钢轨动态求解区中间位置。凹坑深度为0.5 mm，横向宽度为10 mm，纵向长度为10 mm。通过修改钢轨表面局部节点坐标，即可完成上述凹坑缺陷的施加。值得注意的是，钢轨表面施加凹坑缺陷后网格节点会发生位移，需要保证凹陷深度小于0.15 mm(该区域钢轨单元特征尺寸为1 mm)的15%，以免影响轮轨瞬态滚动接触行为^[1]。施加凹坑缺陷后的钢轨表面如图 8所示。

当轨面含有凹坑缺陷时，车轮以55 km/h的平动速度滚过动态求解区，5个代表性时刻(依据图 7选取)所对应的轮轨滚动接触黏滑分布图如图 9所示。

通过对比图 9a和9b与图 7a和7b可知，在车轮滚过凹坑缺陷之前，轮轨接触斑区域黏滑分布与轨面光滑工况下的计算结果相似，表现为接触斑前缘为黏着区域，接触斑后缘为滑移区域；在车轮滚过凹坑缺陷时，轮轨接触斑形状如图 9c所示，表现为环绕着凹坑的2个接触斑区域均为滑移区域；在车轮滚过凹坑缺陷之后，轮轨接触斑区域面积减小，而且接触斑区域几乎表现为滑移状态(尽管图 9d中极少数节点表现为黏着状态，但是绝大多数节点仍表现为滑移状态)，如图 9d和9e所示。综上所述，钢轨表面凹坑缺陷的存在改变了轮轨滚动接触黏滑状态，使轮轨切向力在宏观上有达到饱和的可能性，进而导致轮轨界面发生滑移。轮轨界面一旦发生滑移，就会诱发轮轨系统失稳，进而引起钢轨表面材料磨耗，以致最终可能形成钢轨波磨。以上所述仅说明钢轨表面有发生初始波磨的趋势，如何确保初始波磨进一步发展并形成最终波磨，还需要对轮轨系统的固有振动特性和钢轨纵向磨耗特征展开研究。

4 钢轨波磨形成机理

本文主要通过复模态理论分析轮轨系统的固有振动特性，并依据Archard磨耗模型分析钢轨纵向磨耗特征，以解释实测线路上钢轨波磨的形成机理。

4.1 复模态分析

复模态理论是一种稳态系统分析方法，能够在频域内预测系统的稳定性。本文轮轨系统的运动方程可以表示为^[38]

$ \boldsymbol{M}_{\mathrm{f}} \ddot{\boldsymbol{x}}+\boldsymbol{C}_{\mathrm{f}} \dot{\boldsymbol{x}}+\boldsymbol{K}_{\mathrm{f}} \boldsymbol{x}=0. $

(9)

其中：M_f、C_f和K_f分别表示质量矩阵、阻尼矩阵和刚度矩阵，由于轮轨界面存在摩擦，其均呈非对称特性；$ \mathit{\boldsymbol{\ddot x}}$、$ \mathit{\boldsymbol{\dot x}}$和x分别表示节点加速度向量、节点速度向量和节点位移向量。式(9)的特征方程可进一步表示为

$ \left(\xi^2 \boldsymbol{M}_{\mathrm{f}}+\xi \boldsymbol{C}_{\mathrm{f}}+\boldsymbol{K}_{\mathrm{f}}\right) \boldsymbol{y}=0 . $

(10)

其中, ξ为特征值，y为特征向量。通常，式(10)具有复特征值。根据式(10)，式(9)的一般解可以表示为

$ \boldsymbol{x}(t)=\sum \boldsymbol{y}_j \exp \left(\xi_j t\right)=\sum \boldsymbol{y}_j \exp \left(\left(a_j+m b_j\right) t\right) . $

(11)

其中：(a_j+mb_j)为第m个特征值，m为正整数，a_j和b_j分别为第j阶特征值实部和虚部；y_j为第j阶特征向量；x(t)为节点位移。a_j是估计系统稳定性的重要参数，并由其确定系统的等效阻尼比ζ_j，

$ \zeta_j=-2 a_j /\left|b_j\right|. $

(12)

当a_j≥0时，ζ_j为负，从而x(t)将会随着t增大而增大；同时，系统将趋于不稳定，这表明ζ_j可以作为一个预测不稳定振动发生趋势的系统参数。

基于复模态理论，对在轨面光滑和轨面缺陷工况下的三维轮轨滚动接触有限元模型进行复模态分析。当车轮滚动至钢轨动态求解区中间位置时，执行复模态计算，然后分别提取a_j≥0的系统振型并绘制图 10和11，其中系统振型均进行了变形缩放处理，缩放因子为200。复模态计算关注的频率为0~1 200.00 Hz，且2种条件下提取的振型数量均为46。

由图 10可知，在轨面光滑工况下，轮轨系统出现了2个不稳定振型——第17和第19阶振型，分别对应不稳定振动频率f为322.52和374.20 Hz，对应特征值实部为3.621和4.836，且2个不稳定振型均表现为钢轨相对于轨道板的垂向弯曲振动。根据式(12)可以获得2个不稳定振型的等效阻尼比分别为-0.022和-0.026，负等效阻尼比会引发轮轨系统不稳定振动，进而诱发钢轨波磨。然而，由于车轮经过动态求解区时，在宏观上始终表现为滚动接触，未发生滑移现象，车轮并不会对钢轨表面造成较大磨耗(根据4.2节Archard磨耗模型，轮轨间如果没有相对滑移速度，就不会产生材料磨耗)，所以，钢轨波磨并不容易在轨面光滑工况下形成。

由图 11可知，在轨面缺陷工况下，轮轨系统同样出现了2个不稳定振型——第17和第19阶振型，不稳定振动频率分别为322.57和374.24 Hz，与轨面光滑工况下轮轨系统的不稳定振动频率相近，对应特征值实部为16.798和19.925，且2个不稳定振型与轨面光滑工况下的不稳定振型相似，也表现为钢轨相对于轨道板的垂向弯曲振动。与轨面光滑工况下不稳定振型的特征值实部相比，轨面缺陷工况下不稳定振型的特征值实部明显增大，增加倍数分别为3.639和3.120，且对应的等效阻尼比也分别降低为-0.104和-0.106，这表明轨面缺陷加剧了轮轨系统的固有不稳定振动特性。可以发现，轨面光滑和轨面缺陷工况下的有限元模型复模态分析结果中的不稳定振动频率均处于实测波磨通过频率范围内，同时，由于存在轨面缺陷，在动态求解区，轮轨界面出现了相对滑移现象，因此综合复模态分析结果和轮轨接触黏滑计算结果可知，钢轨波磨的形成机理为轨面缺陷激励引发的轮轨系统的固有不稳定振动，且该不稳定振动表现为钢轨相对于轨道板的垂向弯曲振动。实际上，轮轨表面不可避免地会出现各种缺陷，而这些缺陷均会对轮轨系统造成不利影响，从而加速缺陷的发展并可能诱发二次病害。

4.2 钢轨纵向磨耗特征

为定量分析钢轨波磨的形成过程，并进一步验证本文提出的钢轨波磨形成机理，对轨面缺陷工况下的钢轨纵向磨耗特征进行了分析。在黏着区，轮轨间没有相对滑移速度，故不会产生磨耗，而在滑移区，滑动摩擦力与滑移速度成正比。基于三维轮轨滚动接触有限元模型，通过计算可得轮轨接触区域的法向接触力和相对滑移速度，然后利用Archard磨耗模型便可计算车轮滚过轨面凹坑缺陷后的钢轨纵向磨耗，具体可表示为：

$ V_{\mathrm{r}}=\sum\limits_{\eta=1}^n \Delta V_{\mathrm{r} \eta}=\sum\limits_{\eta=1}^n \frac{k_{\mathrm{r} \eta} \boldsymbol{F}_{z \eta} \boldsymbol{v}_{\mathrm{rw} \eta} \Delta t}{H}=\sum\limits_{\eta=1}^n \frac{k_{\mathrm{r} \eta} \boldsymbol{F}_{z \eta} \boldsymbol{s}_{\mathrm{rw} \eta}}{H}, $

(13)

$ k_{\mathrm{r} \eta}=\left\{\begin{array}{l} \left(31+54 v_{\mathrm{rw} \eta}\right) \times 10^{-6}, \\ \ \ \ \ \ \ \ v_{\mathrm{rw} \eta}<0.2 ; \\ \left(70-10 \frac{\left|v_{\mathrm{rw} \eta}-0.45\right|}{0.25}\right) \times 10^{-6}, \\ \ \ \ \ \ \ 0.2 \leqslant v_{\mathrm{rw} \eta} \leqslant 0.7 ; \\ {\left[40-30\left(v_{\mathrm{rw} \eta}-0.7\right)\right] \times 10^{-6}, } \\ \quad 0.7<v_{\mathrm{rw} \eta} \leqslant 1.0 . \end{array}\right. $

(14)

其中：V_r为钢轨表面材料磨耗体积，m³；η为接触斑内单元序号，η=1，2，…，n；k_rη为无量纲磨耗系数，与轮轨接触压应力和相对滑移速度有关^[39]；F_zη为轮轨法向接触力，kN；z代表轮轨接触法向方向；v_rwη为轮轨接触斑内节点的相对滑移速度，m/s；Δt为磨耗计算时间步长，s；H为钢轨材料的Vickers硬度，GPa；s_rwη为轮轨接触斑内节点的相对滑移量，m。

利用考虑钢轨表面凹坑缺陷激励的三维轮轨滚动接触有限元模型，计算车轮以55 km/h的平动速度滚过凹坑缺陷后的轮轨法向接触力和接触斑内节点的相对滑移量等参数，并基于Archard磨耗模型计算了钢轨纵向磨耗特征。车轮单次通过钢轨表面凹坑缺陷过程所产生的磨耗量如图 12所示。

由图 12分析可知，当车轮经过钢轨表面凹坑缺陷时，会激发瞬态波动磨耗，从而使钢轨表面形成波浪状磨耗，并且波动磨耗随着与凹坑的距离增加而减小，同时，这也反映在车轮经过钢轨表面凹坑缺陷时轮轨系统处于不稳定振动状态，并且不稳定振动随着车轮运行距离增加而减弱。此外，由于钢轨纵向表面不可避免地会连续存在各种尺度的缺陷，因此，在一些实际线路上，发生波动磨耗的范围可能更大。为进一步分析钢轨纵向磨耗的振动特性，对图 12所示的钢轨磨耗量进行了波长谱变换，得到的磨耗量波长谱如图 13所示。

根据图 13可知，钢轨表面纵向磨耗的特征波长为40~50 mm，这与实际线路上的波磨波长测试情况相符，从而在分析钢轨波磨形成过程的同时，进一步验证了本文所提出的钢轨波磨形成机理。

本文从轮轨微观接触和宏观振动2个尺度研究了钢轨波磨的形成机理，研究维度更全面，这区别于已有的仅关注某一尺度的波磨理论研究^{[16, 18-21, 40-41]}。

对于本文采用的隐式+隐式的有限元法，需要在一个处理器为Inter (R) Xeon (R) W-2245 CPU @ 3.90 GHz，机带RAM为64.0 GB的工作站上进行计算，轨面光滑和轨面缺陷工况下的一个完整的轮轨滚动接触过程所需时间分别约为4.6和5.1 d。

5 结论

本文尝试从轮轨瞬态接触黏滑振动出发，探究轨道线路上钢轨波磨的形成机理。首先，根据现场波磨情况，利用有限元软件ABAQUS建立了三维轮轨滚动接触有限元模型并论证了其有效性；然后，分析了车轮运行过程中的接触黏滑状态，并分别讨论了有无轨面缺陷对接触黏滑状态的影响，进而研究了轮轨接触黏滑特性与钢轨波磨形成的关系；最后，依据复模态理论和Archard磨耗模型分析了轮轨系统的固有振动特性和钢轨纵向磨耗特征，从而解释钢轨波磨的形成机理，得出如下结论：

1) 当车轮滚过光滑钢轨时，轮轨接触斑区域均表现为接触斑前缘为黏着区域，中后缘为滑移区域，这与采用稳态动力学计算的结果相近，表明建立的有限元模型有效。同时，轮轨接触未达到全滑移状态，即车轮一直处于稳定滚动状态，说明在轨面光滑工况下轮轨系统处于稳定状态，因此不易形成钢轨波磨。

2) 当车轮滚过轨面凹坑缺陷时，轮轨接触斑区域表现为环绕凹坑的2个滑移区域；当车轮滚过凹坑缺陷之后，轮轨接触斑区域面积减小，且接触斑区域几乎处于滑移状态。凹坑缺陷改变了轮轨滚动接触黏滑状态，促使轮轨界面发生滑移，进而诱发轮轨系统失稳，并引起钢轨表面材料磨耗，以致最终可能形成钢轨波磨。

3) 复模态分析表明轨面缺陷加剧了轮轨系统的固有不稳定振动特性，且不稳定振动频率均处于实测波磨通过频率范围内。结合轮轨接触黏滑计算结果可知，钢轨波磨的形成机理为轨面缺陷激励引发的轮轨系统的固有不稳定振动，且该不稳定振动表现为钢轨相对于轨道板的垂向弯曲振动。

4) 当车轮经过轨面凹坑缺陷时，会激发瞬态波动磨耗，从而使钢轨表面形成波浪状磨耗。钢轨表面纵向磨耗的特征波长为40~50 mm，这与实际线路上的波磨波长测试情况相符，从而在分析钢轨波磨形成过程的同时，进一步验证了钢轨波磨的形成机理。