高速铁路轮轨垂向力与轨道高低不平顺关联特性

牛留斌

（中国铁道科学研究院集团有限公司 基础设施检测研究所，北京 100081）

良好的轨道平顺性是车辆安全运行的前提，不良的轨道状态会导致轮轨间产生异常作用力，较大的轮轨力引起车辆-轨道系统强烈振动冲击和疲劳损伤，降低车辆运行安全。作为辅助评价轨道服役状态的重要参数，国内外铁路机构采用测力轮对技术实时测量轮轨间的作用力［1-3］，评估车辆运行安全风险状态，为轨道线路的养护维修提供参考依据。

轮轨力是轮轨耦合作用最直接的结果，其大小受到车辆、轨道、轮轨接触状态等多种因素的影响，如轨道高低不平顺、车轮型面磨耗状态［4］、踏面不圆顺度［5］、道岔区结构［6］、轨道支撑刚度［7］、道床刚度突变［8］和轮轨之间的摩擦系数［9］等。作为车辆-轨道系统的主要激扰源，轨道不平顺是引起轮轨力变化的重要因素。以轨道不平顺作为轨道-车辆系统的输入量，研究轨道不平顺引起的轮轨力、行车安全、车辆振动状态之间关系的课题相对较多，如Choi 等［10］采用VAMPIR 软件，仿真研究直线段轨道高低、水平、轨向等轨道不平顺引起的轮轨力波形变化特征和运行安全的影响；Steenbergen［11］采用有限单元方法和Timoshenko梁轨道模型，研究焊接接头处轨道短波不平顺与轮轨力响应之间的时频分布特征；Seo 等［12］采用商业软件建立三维有限元模型，仿真计算轨顶塌陷深度、轨面裂纹长度等轨道状态条件下的轮轨力响应与发展趋势，并讨论牵引力与制动力对轨面裂纹发展速度的影响；Pombo 和Naeimi［13-14］等采用多体动力学软件与轮轨耦合动力学理论，仿真研究小曲线超高、水平不平顺等轨道因素对车辆运行安全的影响规律及相应轮轨力变化特点；Tanaka 等［15］采用系统辨识与数据建模方法，利用传递函数和10 m弦测轨道高低不平顺数据反演轮轨垂向力和轴箱加速度；Gullers 等［16］采用实测轮轨力数据对轨道几何不平顺和轨道刚度变化地段进行状态评价和现场复核，总结轨道状态与轮轨力波形之间的对应关系。轨道高低不平顺与轮轨垂向力之间关联特性一直是轮轨关系的研究重点，波长和幅值是轨道高低不平顺的主要特征参数，上述文献侧重研究了轨道不平顺幅值、变化梯度等与轮轨力之间的对应关系，而波长与车辆动力性能方面的研究侧重于车辆敏感波长［17-18］分布特征和轨道长波不平顺管理限值建议，对轮轨耦合作用下轨道高低不平顺波长与轮轨垂向力关联关系的研究相对较少，二者关系的深化研究有助于掌握波长特性激励的轮轨动力学行为，完善轨道不平顺波长管理体系。

本文以我国高铁线路和典型服役车辆参数建立三维轮轨接触有限元模型，分析高速条件下轨道高低不平顺幅值、波长与轮轨垂向力在时频域上的关联关系，拟合附加动态轮轨垂向力与轨道高低不平顺波长之间的函数表达式。结合以上研究，提出了利用轨道高低不平顺反演附加动态轮轨垂向力算法及算例，为我国线路精准养护维修提供支撑依据。

1 有限元模型构建

1.1 模型构建

车辆在钢轨上运行时，轮轨间作用力可用静轮重和附加轮轨垂向力［19］表述为

式中：F为轮轨垂向力，kN；F0为静轮重，kN；Fdyn为轨道高低不平顺引起的附加轮轨垂向力，kN。

钢轨在轮轨垂向力F作用下产生弹性挠曲变形，轨道静态高低不平顺发展成轨道动态高低不平顺，轨道变形示意图如图1所示。

图1 轨道在轮轨垂向力作用下变形示意图

轨道高低不平顺激励下轮轨三维有限元模型结构简图如图2 所示。图中，坐标系原点位于前轮轨接触O点位置，轮对运行方向为z轴方向，车辆一系悬挂及以上部件简化为刚度K1、阻尼C1的弹簧单元和质量块m，用等间隔的弹簧单元（刚度K2、阻尼C2）模拟扣件系统的支承作用，轴箱等轮对附属部件的质量均匀施加在轮对上，轨道板、CA砂浆层、混凝土道床和路基等统一简化为基础支承单元。

图2 轮轨接触有限元模型结构简图

轮轨三维有限元模型如图3 所示。模型中，车轮和钢轨均采用其真实几何尺寸，车轮踏面为LMA型，钢轨廓形为CN60型，轨道底部设置1∶40轨底坡；轨枕间距为650 mm，由弹簧阻尼单元等效，钢轨长度为41.6 m，包含65组等效弹簧单元，前后轮对轴距为2.5 m。为模拟轮轨接触的真实应变状态，采用双线性弹塑性本构关系模型模拟轮轨材料的弹塑变形性能。在轮轨接触面的法线方向采用“面-面”硬接触算法计算轮轨法向接触力，切线方向采用干燥条件下的轮轨相互作用状态，且库伦摩擦系数设为0.4。轨道端部施加面对称约束，为减少轨道端部约束对计算结果的影响，轮对初始位置与钢轨端部的距离不小于3 m，轨道高低不平顺设置在图2 所示的钢轨中部ab段。为缩减有限元模型计算时间，单元网格尺寸与单元所在区域的关注程度有关，靠近轮轨接触区域的网格尺寸较小，如在轮轨接触面上实体单元网格的最小尺寸为2 mm，路基等远离轮轨接触区域的单元网格尺寸较大。

图3 轮轨三维有限元网格模型

为减少ABAQUS 显式模块在初始时刻动态平衡求解所需的计算时间，先采用隐式模块计算重力载荷下模型的静平衡状态，并将该平衡状态作为显式模块计算初始时刻的边界条件。在初始时刻轮对上施加加速度载荷，轮对到达图2 所示钢轨a点时有限元模型处于动态平衡状态。为减少轮轨运行过程中轮轨间摩擦力对车辆运行速度的影响，在车轮上施加牵引扭矩M以使轮对匀速运行在轨道上。牵引扭矩M由轮对转动定理公式［19］得到，为

式中：Fz为轮轨作用力F沿车辆运行方向即z轴方向上的分量，kN；R为轮对滚动圆半径，430 mm；J为轮对的转动惯量；ε为轮对角加速度，车轮匀速行驶时ε=0。

在不考虑轮对横向运动的前提下Fz可由下式计算。

式中：μ为牵引系数，下文取0.03。

ABAQUS 显式模块采用中心差分法显式积分求解轮轨实时运动状态和接触力，每个时间增量步结束时刻的状态完全由该增量步开始时刻的位移、速度、加速度计算得到［20］。最大积分时间步长Δt由有限元模型最高固有频率ωmax决定，满足下式。

其中，

式中：ξ为轮轨系统的临界阻尼比；Le为模型单元长度，mm；cd为模型材料特性决定的波速，m·s-1；E和ρ分别为材料的弹性模量和密度。

由式（4）可知，图3 所示轮轨有限元模型中网格的最小尺寸为2 mm 时，最大积分时间步长Δt的数量级为10-7s。

1.2 数据处理步骤

在图3 所示有限元模型中的钢轨上施加轨道高低不平顺，以研究其引起的轮轨垂向力变化特性。施加的轨道高低不平顺来源于中国高速铁路无砟轨道高低不平顺谱样本［21］，为

式中：f为空间角频率，m-1；Sv(f)为轨道高低不平顺谱，mm2·m；C和B为分区段拟合系数，rad·m-1，具体参数取值见表1。

表1 高速铁路无砟轨道高低不平顺谱拟合系数

通过修改钢轨横断面上全部网格节点坐标的方式，在图2 中ab 段轨道模型上施加式（5）生成的轨道高低不平顺样本。施加轨道高低不平顺时模型网格修正示意图如图4所示。

在图4 中钢轨网格节点G处施加幅值为dzG的轨道高低不平顺时，将原节点G处的坐标（xA，yA，zA）修改为新节点G′处的坐标（XG′，YG′，zG′），换算公式如下。

图4 施加轨道高低不平顺时模型网格修正示意图

式中：θ为1∶40 轨道坡角度；h为钢轨高度，CN60 型钢轨h值为176 mm；za和zb分别为图2 中钢轨a点和b点截面位置处纵坐标。

综上，不同工况下轨道高低不平顺仿真计算及数据处理步骤如下。

（1）建立轮轨三维有限元模型，在钢轨ab段施加轨道高低不平顺。

（2）隐式计算静载荷作用下轮轨位移场，作为显式计算时轮轨模型初始边界状态。

（3）施加运行速度和牵引扭矩，显式计算轮对经过轨道高低不平顺时的动力学响应。

（4）按照4 mm 间距依次读取钢轨ab段单元节点上轮轨垂向力时间历程数据，输出频率为10 kHz。

（5）轨道高低不平顺区段上各节点轮轨垂向力时间历程数据低通滤波。

（6）修改钢轨ab段单元节点坐标，以模拟施加轨道高低不平顺，仿真计算并输出新工况下轮轨垂向力。

（7）汇总各个工况下轮轨垂向力结果，分析高低不平顺波长、幅值与轮轨垂向力变化规律，趋势拟合等。

2 模型验证及数据处理

除特别标出外，图3 所示三维有限元模型的仿真参数主要源于我国典型高铁线路及服役车辆，取值见表2。

表2 轮轨有限元模型主要参数

车辆运行速度为300 km·h-1且未施加轨道高低不平顺时的轮轨垂向力波形如图5 所示。由图5可知：轨枕不连续支承作用引起轮轨垂向力在静轮重70.3 kN 附近周期性波动，波动幅度为±0.5 kN，波长等于轨枕间距。

图5 未施加轨道高低不平顺时轮轨垂向力波形图

利用式（6）将余弦型轨道高低不平顺施加到图4有限元模型中钢轨ab段，得到轨道高低不平顺波长为1.5 m、幅值为1 mm 时轮轨垂向力波形，并与文献［22］在相同参数下的仿真结果进行对比，如图6所示。

图6 轮轨垂向力仿真波形及对比

由图6 可知：未施加轨道高低不平顺区段内的轮轨垂向力波动，是由不连续的轨枕支承弹簧引起的；本文轮轨垂向力仿真波形与文献［22］数值结果吻合，变化趋势一致，验证了图3 所示三维有限元模型仿真结果的正确性。

我国常规轨道高低不平顺波长范围为1.5～120.0 m［23］，但限于仿真运算的速度和规模，下文计算了1.5～12.0 m 波长范围内轨道高低不平顺引起的轮轨垂向力特征，波长大于12 m 的轨道高低不平顺引起的轮轨垂向力由上述仿真结果拟合推理得到。轮轨垂向力数据在分析之前，先对数据进行低通滤波处理。图6 中轮轨垂向力大值随低通滤波截止频率变化的曲线如图7 所示。由图7 可知：当截止频率小于1 500 Hz时，轮轨垂向力大值随截止频率的增加而快速增大；当截止频率在1 500～2 500 Hz频段时，轮轨垂向力大值随截止频率增速变缓；当截止频率在2 500～5 000 Hz 频段时，轮轨垂向力大值趋于稳定。因此，在不同高低不平顺激励的轮轨垂向力数据分析之前，先对各工况的仿真结果进行2 500 Hz低通滤波处理。

图7 轮轨垂向力大值随截止频率变化曲线

3 数值结果

按式（6）生成的中国高速铁路无砟轨道高低不平顺样本数据如图8所示。在仿真计算时将轨道高低不平顺的幅值进行放大而波长保持不变，以研究不同幅值轨道不平顺对轮轨垂向力的变化特效。幅值放大系数（Amplitude Multiple Factor，AMF）分别为2，4，5 和8 时轮轨垂向力仿真结果如图9所示。

图8 中国高速铁路无砟轨道高低不平顺样本谱

由图9 可知，虽然放大系数AMF 由1 增至8，但轮轨垂向力波形、相位等变化趋势相似，只是轮轨垂向力局部最值与静轮重之间的差值增大，即轮轨垂向力波动范围变大。

图9 不同AMF条件下轮轨垂向力波形

定义幅值增大倍数（Amplitude Increase Factor，AIF）为不同放大系数时仿真结果与AMF=1时的差值，不同AIF 条件下轮轨垂向力波形如图10 所示。图中：虚线、实线与波形相交位置的纵坐标标值分别为轮轨垂向力的最大值、最小值。

图10 不同AIF条件下轮轨垂向力波形

图10 中，不同列（从左至右、先虚线后实线顺序编号）的轮轨力分别见表3和表4。

表3 轮轨垂向力局部最大值

图9 中在AMF 为8 时，部分轮轨垂向力为0 kN，表明此时轮轨已脱离接触，所以表3 中第1，3，6，8—11 列与表4 中各列轮轨垂向力数值之间的比例约为1∶3∶4∶7，与AIF值相近。

表4 轮轨垂向力局部最小值

由以上可知，轨道高低不平顺幅值的变化影响附加轮轨垂向力的大小，但不影响其波形相位。所以在研究轮轨垂向力与轨道高低不平顺波长的关系时，可先将其幅值设为1 mm。轨道高低不平顺波长λ=1.5 m、幅值为1 mm 时的轮轨垂向力仿真结果如图11 所示。由图11 可知：轮轨垂向力的时域波动周期与轨道高低不平顺的激励周期相等；轮轨垂向力在静轮重F0附近呈现振幅为α的周期性变化，与轨道高低不平顺波形之间存在相位差（延迟相位）β。

图11 波长1.5 m轨道高低不平顺条件下轮轨垂向力仿真结果

图11 中轮轨垂向力波动成分即为轨道高低不平顺激励的附加轮轨垂向力Fdyn，其频谱曲线如图12 所示。图中：圆圈处为轨道高低不平顺引起的Fdyn谐波成分，五角星处为轨枕引起的谐波成分。

图12 轮轨垂向力频谱曲线

由图12 可知：与Fdyn频谱曲线的主频相比，2阶及以上谐波成分能量均较小，在Fdyn的总能量中所占比重不足3%；主频对应的空间频率为0.667 m-1，对应的空间波长为1.5 m，即Fdyn主频是由轨道高低不平顺波长决定的。

因此，幅值为A、波长为λ的轨道高低不平顺引起的Fdyn，λ主要特征参数为α，β和λ，可记为

其中，

式中：z为图2 所示车辆沿着钢轨运行过程中轮轨接触位置的纵坐标；α和β为幅值1 mm、波长λ的轨道高低不平顺引起附加轮轨垂向力的振幅和延迟相位；ϖλ为轨道高低不平顺的激励角频率；ν为车辆运行速度。

在轮轨三维有限元模型中施加不同波长且幅值为1 mm 的轨道高低不平顺，得到不同工况下附加轮轨垂向力波形的振幅α、延迟相位β与波长λ之间关系曲线如图13所示。

图13 不同波长工况下附加轮轨垂向力波动幅值与延迟相位曲线

由图13 可知：波长λ从1.5 m 增加到10.0 m时，附加轮轨垂向力的波动振幅α由39.5 kN 减少到1.5 kN，而延迟相位β则 由0.087π 减少到0.015π，且其分布特征与钢轨的垂向振动特性及激励的频率有关。

图3 所示有限元模型中钢轨垂向位移导纳曲线如图14所示。曲线中第1个共振峰值对应钢轨共振频率frail，其值可用下式估算［2，24］。

式中：k为单位长度钢轨的轨下刚度，N·m-2；mr为单位长度的钢轨质量，kg·m-1。

波长λ为1.5～120.0 m 的轨道高低不平顺在车速为300 km·h-1时激励频率约为0.7～56.0 Hz，该频段对应图13 中黄色区域远离钢轨的共振频率frail，且该频段的钢轨垂向导纳曲线未出现大的共振峰线。波长λ较大时，轨道高低不平顺激励的附加轮轨垂向力频率较小，由图14 所示钢轨垂向位移导纳曲线可知钢轨动态效应也较小，附加轮轨垂向力与轨道高低不平顺几乎同时发生，二者延迟相位β较小；波长λ较小时，轨道高低不平顺引起附加轮轨垂向力频率较高，钢轨的动态效应增强，二者之间延迟相位β增大。

图14 轮轨有限元模型中钢轨垂向位移导纳曲线

利用非线性最小二乘方法［25］对图13中附加轮轨垂向力波动的振幅α及延迟相位β进行有理式曲线拟合，得

式中：p1—p4和q1—q3为待定拟合参数。

在拟合优度大于99.5%时，α，β与波长λ之间的拟合式及参数取值如图15 所示。图中：λ＞12 m 的轨道高低不平顺激励的附加轮轨垂向力的α和β值依据上述仿真结果进行拟合外推，即图中拟合曲线的右侧实线段。

轨道高低不平顺引起的激励频率还与车辆运行速度有关，车速提高至350 km·h-1条件下轨道高低不平顺引起附加轮轨垂向力的振幅α及延迟相位β与波长λ之间的拟合方程如图16所示。

对比图15 与图16 可知：车速350 km·h-1时较300 km·h-1时，振幅α的增加幅度明显大于延迟相位β的增加幅度，如波长为1.5 m 时振幅α由300 km·h-1时的39.5 kN 增大到350 km·h-1时的54.8 kN，增加了38.7%，而延迟相位β则由0.087π 增大到0.095π，增加了9.2%；随着波长λ增加，振幅α和延迟相位β分别趋近于1λ和0.04π。

图15 车速300 km·h-1时轨道高低不平顺波长与附加轮轨垂向力振幅、延迟相位拟合曲线

图16 车速350 km·h-1时轨道高低不平顺波长与附加轮轨垂向力波动振幅、延迟相位拟合曲线

将式（9）和式（10）代入式（8），得幅值A、波长λ轨道高低不平顺激励下的附加轮轨垂向力Fdyn，λ为

定义轨道变化率η为轨道高低不平顺幅值A与波长λ比值，代入式（13）有

其中，

由式（14）及图14 拟合结果，得到车速300 km·h-1时附加轮轨垂向力大值Fdyn，λmax与轨道变化率η之间的斜率k斜率随轨道高低不平顺波长λ变化曲线如图17所示。

图17 附加轮轨垂向力大值和轨道变化率之间斜率随波长变化曲线

由图17 可知：当轨道高低不平顺波长λ大于30 m 时，斜 率k斜率中Fdyn，λmax与η之间的斜率逐步接近式（9）中参数p1，此时Fdyn，λmax为

由式（15）及图17 可知，在轨道高低不平顺波长λ大于30 m 时，附加轮轨垂向力最大值与轨道变化率η 近似呈现线性关系。在表2 所列车辆轨道参数条件下，车速300 km·h-1时p1为13.14；车速350 km·h-1时p1为19.19。

高铁线路上的轨道高低不平顺包含众多波长成分，此时附加轮轨垂向力Fdyn为单个波长成分引起的附加轮轨垂向力的和，可表示为

式中：i为轨道不平顺波长成分数量；λi为第i个轨道高低不平顺成分的波长；Fdyn，λi为第i个轨道高低不平顺成分引起的附加轮轨垂向力。

由式（16）及傅里叶变换［26］可由轨道高低不平顺反演得到附加轮轨垂向力，步骤如下。

（1）对实测轨道高低不平顺s(λ)数据进行快速傅里叶变换，得到空间频域上的傅里叶变换对s(ω)，即

式中：F 为傅里叶变换符号；ω为波长的倒数，即空间频率。

（2）将各频段轨道高低不平顺傅里叶变换对s(ω)的幅值乘以附加轮轨垂向力的振幅α、减去相位延迟量β，得到修正后的轨道高低不平顺傅里叶变换对S(ω)，即

（3）对修正后的轨道高低不平顺S(ω)进行傅里叶逆变换，得到附加轮轨垂向力F(λ)，即

按照上述步骤，采用某高铁线路车速300 km·h-1区段实测轨道高低不平顺反演得到了附加轮轨垂向力波形，并与该区段实测轮轨垂向力数据进行对比，结果如图18 所示。其中，实测轨道高低不平顺的波长范围为1.5～120.0 m，按照该波长范围和车辆运行速度对实测轮轨垂向力数据进行带通滤波处理。

由图18 可知，采用实测轨道高低不平顺和上述步骤反演得到的附加轮轨垂向力与其实数据吻合良好，变化趋势一致，验证了上文研究及拟合式即式（10）—式（12）的正确性。

图18 利用轨道高低不平顺反演附加轮轨垂向力与实测轮轨垂向力对比

4 结 论

（1）为研究高速条件下轨道高低不平顺激励的轮轨垂向力特性，采用ABAQUS 软件及轮轨真实尺寸建立轮轨接触有限元模型，通过与对比相关文献仿真结果，验证所建模型的准确性。在钢轨模型中输入中国高速铁路无砟轨道高低不平顺谱样本数据，研究轨道高低不平顺幅值对轮轨垂向力的影响特性。结果表明：在轮轨不脱离接触的前提下，附加轮轨垂向力波动变化量与轨道高低不平顺的幅值成正比。

（2）仿真计算了车速300 km·h-1时轨道高低不平顺波长为1.5～12 m、幅值为1 mm 轮轨动力学响应。轨道高低不平顺引起的一阶谐波占附加轮轨垂向力的主要成分，其余谐波成分所占比重小于3%；附加轮轨垂向力的波动振幅随着轨道高低不平顺波长增加而非线性减少，二者的对应关系可用有理式方程进行拟合，轮轨力波形与轨道高低不平顺之间相位延迟量随激振频率的增大而增加，整体处于0.02π～0.10π之间。

（3）文中给出了车速300 和350 km·h-1时轨道高低不平顺幅值、波长与附加轮轨垂向力之间的函数表达式。附加轮轨垂向力最大值与轨道高低不平顺幅值、波长之间不存在明显的线性关系，但轨道高低不平顺波长大于30 m 时，附加轮轨垂向力最大值与轨道高低不平顺的变化率近似线性关系，二者的比例系数与车辆、轨道参数及车速等因素有关，车速为300 和350 km·h-1时比例系数分别为13.14 和19.19。由高铁线路上实测轨道高低不平顺反演的轮轨垂向力结果与实测轮轨垂向力比对验证了拟合方程表达式的正确性。