随着“一带一路”倡议的开展以及世界铁路通车里程的日趋增加，铁路正改变着人们的出行方式，而轮轨相互作用产生的振动和噪声越发影响着乘车的舒适性与安全性. 为改善轮轨相互作用，多种减振措施被提出并投入使用，钢轨扣件具有提供系统弹性、连接钢轨与轨枕以及减缓振动、延缓轨道系统损伤破坏等作用^[1]. 高分子橡胶材料制成的弹性垫板在钢轨扣件中起到至关重要的作用，而弹性垫板在受到动态应力及应变时会产生一定记忆性变形及能量耗散，致使在长期服役情况下，钢轨扣件橡胶垫板的刚度/阻尼具有复杂的非线性力学特征，其刚度具有动态荷载频率与幅值相关性^[2-4].

早期的车辆-轨道耦合动力学模型通常利用Maxwell流体模型将扣件简化为线性弹簧和黏性阻尼串联的形式，或利用Kelvin-Voigt固体模型把扣件考虑为线性弹簧和黏性阻尼并联的形式，2种简化模型改变了阻尼因子及刚度弹簧进而影响动态轮轨相互作用，而后发展到广义Maxwell模型和广义Kelvin-Voigt模型，但计算参数过于复杂不适用于工况计算^[4-7]. 当前Kelvin模型因结构简单、计算高效而在车辆-轨道耦合动力学模型中被广泛使用，但其未考虑频率及幅值变化引起的胶垫非线性力学特征，导致计算精度不足，与实际运营情况不符.

由此可见，如何正确表征扣件橡胶垫板的本构关系及复杂的非线性变化是精确仿真结果的重要环节. 对此许多学者进行了深入研究：Gil-negrete等^[3-4]验证Kelvin-Voigt模型属性并建立模拟橡胶垫板的五参数模型；Berg等^[8]利用微积分的方法描述橡胶材料受荷非线性变化并建立一维非线性模型；Zhu等^[9]建立了车辆-轨道耦合动力学模型并考虑了轨道板的频率幅值依赖性；张大伟等^[10]建立车辆轴箱非线性橡胶垫板并分析其动态相互作用.

为了更加精确地评估动态荷载作用下弹性橡胶垫板非线性行为及动态黏弹性力学特征，进一步提高三维轮轨滚动接触瞬态动力学有限元模型的准确性，本文以我国高速铁路无砟轨道扣件WJ-8弹性橡胶垫板为实验研究对象，通过静刚度实验反演非线性动力学相关参数，并运用ANSYS/LS-DYNA大型有限元软件建立充分考虑弹性橡胶垫板非线性变化及动态黏弹性特征的钢轨-扣件波磨模型，分析在波磨作用下，橡胶垫板非线性特征对轮轨接触高频响应的影响.

1.   扣件橡胶垫板黏弹性力学性能

1.1   扣件橡胶垫板力学模型

本文建立的动力学模型充分考虑了扣件橡胶垫板复杂非线性变化与动态黏弹性力学性能，由传统的弹性因子、并联分数阶导数黏弹性因子以及非线性摩擦因子3部分组成^[1]，相较传统的瞬态动力学模型，其更加真实地反映了轮轨相互作用和扣件橡胶垫板的频率相关性及幅值相关性，所用模型如图1所示. 对试件扣件施加的荷载为

图  1  非线性扣件橡胶垫板力学模型

Figure  1.  Mechanical model of rubber pad of nonlinear fastener

下载: 全尺寸图片 幻灯片

式中：${F_{\rm{e}}}\left( t \right)$为传统的弹性力，t为时间；${F_{\rm{v}}}\left( t \right)$为与振动频率具有相关性的分数阶导数黏弹性力；${F_{\rm{f}}}\left( t \right)$为与振动幅值具有相关性的非线性摩擦力.

传统的弹性力由弹簧刚度K_e与位移x乘积得到，如式（2）.

橡胶垫板的黏弹性行为受到长期服役条件的影响，具有一定的记忆性，随服役时间而产生疲劳变形和应力松弛，对Bagley-Torvik^[11]分数阶微分方程进行推广，使方程能够较好地反映出橡胶垫板黏弹性行为的时间记忆性，分数阶导数黏弹性力的时域离散公式为

式中：${F_{\rm{v}}}\left( {{t_N}} \right)$为橡胶垫板分数阶导数黏弹性力；$t_{{N}}$为对时间离散化处理后的第N个时间段长度；${C_{\alpha } }$为阻尼系数；$_0^GD_t^\alpha$ 为微分算子，左上角字母G代表算子类型，右上角字母$\alpha$代表微分阶次，左下角0代表积分起点，右下角$t$代表积分终点，D为扣件橡胶垫板的正向位移值；$\Delta t$为微小的时间段；${x_{{N}}}$为第N个时刻的橡胶垫板位移值；$A_{j + 1}$为第j+1个时刻的Grunwald系数；${x_{N - j}}$为第N−j个时刻的橡胶垫板位移值.

根据橡胶增强理论，当前钢轨扣件橡胶材料多为合成橡胶，继而在橡胶材料受荷时产生复杂的摩擦接触和摩擦力，摩擦力具有非线性特征，其值与振动幅值具有一定关系，因此，在力学模型中并联一个非线性摩擦因子，非线性摩擦力与位移的关系为

式中：${F_{{\rm{fs}}}}$为位移达到状态参考点${x_{\rm{s}}}$时的摩擦力；$\mu = {F}_{{\rm{fs}}}/{F}_{{\rm{f}}\mathrm{max}}\in \left[{-1}，\text{1}\right]$，${F_{{\rm{f}}\max }}$为最大摩擦力；$x$为线弹性变形对应位移值，此时代表摩擦力下对应的位移值，${x_2}$为每个滞回环摩擦力对应位移的中间值，$\left( {{x_{\rm{s}}},{F_{{\rm{fs}}}}} \right)$为滞回环上每个分支的状态参考点.

1.2   力学模型参数确定

本文中模型参数取值参考文献[1]，文献[1]中试验是利用万能试验机和高速铁路无砟轨道扣件WJ-8完成的，试验台上自下而上依次平整地摆放支承钢板、弹性橡胶垫板、加载钢板，为了使试验垫板受力更加真实，在加载钢板上部摆放0.5 m长的CHN60型短钢轨，试验中，垂直预加载100 kN，卸载后等待1 min，再次加载100 kN，卸载后等待1 min，然后将位移指示器调零进行正式试验，以2～3 kN/s的速度对试验短钢轨垂直加载，根据式（5）求解试验扣件静刚度.

试验选取3组WJ-8扣件，每组扣件进行两次加载试验，对试验结果进行数值拟合，寻找规律以反推钢轨扣件橡胶垫板力学模型的试验参数，试验所得橡胶垫板静刚度与加载值关系如图2^[1]所示.

图  2  试验测得扣件静刚度与荷载拟合曲线

Figure  2.  Fitting curve between static stiffness and load of fasteners measured by test

下载: 全尺寸图片 幻灯片

根据试验结果反推力学模型待定参数，计算中，首先输入初始条件扣件弹条扣压力18 kN，然后从0.01 Hz的频率施加准静态荷载0～70 kN，0.5 m短钢轨质量为30.32 kg，计算得到${C_{\alpha }}{\text{ = 6.2}} \times {\text{10}}^{\text{6}}\;{{\rm{Ns}}^\alpha }/{\rm{m}}$，$\alpha {\text{ = 0}}{\text{.5}}$，F_{f max}=1 200 N，${x_2}{\text{ = 2}} \times {\text{1}}{{\text{0}}^{{{ - 4}}}}\;{\rm{m}}$，${K_{\rm{e}}} = 1 \times {10^6}\; {\rm{N/m}}$，试验结果与理论具有一致性，橡胶垫板受荷载后其刚度随荷载增加而降低，静刚度变化具有明显的非线性特征.

由此可知，以往的三维瞬态滚动接触研究中，未考虑钢轨扣件橡胶垫板的刚度非线性特性，导致计算结果精度不足，本文根据以上试验数据支撑，将橡胶垫板的非线性变化及动态黏弹性力学特征引入三维瞬态滚动接触有限元模型中，以提升计算结果的精确性.

2.   三维瞬态滚动接触有限元模型

2.1   模型概述

应用ANSYS/LS-DYNA有限元软件参考文献[12]建立钢轨波磨三维瞬态有限元模型，考虑真实车轮与钢轨的几何形状、轮轨高频响应以及柔性轮轨，基于三维轮轨瞬态滚动接触显式有限元模型在时域范围内的真实受力特征，反映车轮沿基本轨单向驶过波磨区段的轮轨相互作用. 根据车辆-轨道系统关于钢轨纵向中心线具有对称性，所建模型为半轮对及单侧钢轨，整个仿真模型包括弹簧上质量、一系悬挂、车轮、钢轨以及轨下扣件系统，如图3所示，车轮选用LMA型踏面，钢轨选用CN60轨，扣件间距0.65 m，为消除初始激扰对求解的影响，在初始位置至求解区设置1.50 m动力松弛区，波磨求解区为0.85 m，为消除边界效应影响，在车轮滚动区域两侧各设置5.00 m边界，模型总长12.35 m.

图  3  三维轮轨瞬态滚动有限元模型

Figure  3.  Finite element model of three-dimensional wheel–rail transient rolling

下载: 全尺寸图片 幻灯片

为精确模拟高速列车在运营条件下驶过波磨区段的过程，对不同行驶速度以及不同波磨深度的工况进行计算，进一步揭示考虑橡胶垫板频变特性后速度及波深对轮轨相互作用的影响，采用8节点六面体实体单元离散真实几何形态的轮轨，采用考虑真实车体及转向架质量的质量单元模拟簧上质量，一系悬挂及轨下扣件系统采用平行的弹簧-黏滞阻尼单元阵列，一系悬挂计算参数利用K-V模型固定刚度及阻尼系数，扣件系统采用非线性模型，刚度系数及阻尼系数取自试验测得的真实非线性刚度及阻尼，计算参数如表1所示.

表  1  模型参数

Table  1.  Model parameters

簧上质量/kg	车辆一系悬挂			簧下质量			车轮及钢轨材料
	刚度/ （MN·m−1）	阻尼/ （kN·s·m−1）		车轮质量/kg	轮下附属部件质量/kg		弹性模量/ GPa	密度/ （kg·m−3）	泊松比	阻尼常数
6000	0.88	4		656	340		210	7790	0.3	0.0001

下载: 导出CSV 

| 显示表格

为尽量减小计算规模又可以精确表征胶垫频变特性对轮轨动态相互作用及接触行为的影响规律，根据Lagrangian描述进行不均匀网格划分，轮轨接触区域精细划分为1 mm网格，远离接触区域的网格逐渐变粗，模型划分单元总数1249244，节点总数1430426. 采用双线性随动强化模型进行轮轨本构关系模拟，充分考虑了车轮-钢轨系统在传力过程中的非线性影响关系. 轮轨间滚动接触采用面-面罚函数接触算法表征，利用ICONT接触单元实常数来描述车轮与钢轨在运行过程中的高频相互作用. 计算步骤忽略车轮的横向位移及自旋运动，即对车轴进行横向及纵向约束，一系悬挂进行相同约束，扣件系统底部及钢轨两端全约束，对车轮施加重力场，利用ANSYS在隐式计算下求得节点位移，将节点位移代入LS-DYNA，拆除车轮及簧上质量的纵向约束，并施加纵向初始速度及转动力矩，基于显式时间积分算法求得车轮滚过波磨求解区的瞬态动力响应.

2.2   波磨施加

利用修改求解区表面节点坐标的方式施加规则余弦函数波磨几何形状，根据调研，当节点坐标修改量远小于表面单元网格尺寸时，节点坐标修改对计算结果影响可忽略. 某高速铁路波磨现场实测数据，该段铁路波磨波深主要集中于120～240 μm，因此以波深120 μm为例，求解区钢轨波磨沿纵向施加示意如图4.

图  4  波磨施加示意

Figure  4.  Corrugation application

下载: 全尺寸图片 幻灯片

3.   模型验证与结果分析

3.1   扣件弹性垫板非线性力学特征对轮轨力影响

3.1.1   时域内轮轨力结果分析

本文基于考虑扣件橡胶垫板非线性力学特征的三维瞬态滚动接触有限元模型，对高速列车驶过真实波磨几何形状钢轨进行模拟仿真，分别对不同运行速度（250、300、350 km/h）及不同波磨波深（120、180、240 μm） 2种工况仿真计算并提取轮轨力、轴箱加速度等信息进行时频域分析.

提取轮轨竖向力时程曲线幅值及变化趋势与文献[12]基本一致，数值不同是设置不同的运行速度和波磨波磨形式引起，并在扣件位置处出现同样的轮轨动态力偏大现象，由此验证考虑非线性扣件橡胶垫板的三维瞬态滚动接触有限元模型的正确性. 为验证模型精度，根据波磨实验测试数据总结不同行车速度对应激振频率谱^[13]，行车速度为250 km/h时，列车行驶至波磨区段产生激振频率为500～1000 Hz，多集中于800 Hz处，与模型计算所得850 Hz相近， 350 km/h时，激振频率为650～1560 Hz，多集中于1050 Hz，与模型计算所得1050 Hz一致，模型在高速运行工况下精度更高，因此，采用本文优化完善的模型进行后续瞬态动力学研究是可靠有效的.

图5展示了运行速度250 km/h时橡胶垫板刚度非线性对不同波磨深度（波深）工况时域内的轮轨竖向力影响，并给出相应线弹性轮轨竖向力时程曲线. 由图5可知：轮轨力随波深的增加而更加剧烈，波深由120 μm增长到180 μm时，轮轨竖向力增长幅度最大可达到18.12%，切向力增幅25.69%；当波深由180 μm增长到240 μm时，轮轨竖向力增幅最大为29.68%，切向力增幅47.12%，显著增加的轮轨相互作用与波磨发展后，车轮行驶至波磨波峰位置的挤压碰撞有关；仿真计算表明当波磨深度达到240 μm时，车轮滚过波磨求解区所受轮轨竖向力最大值为144.553 kN，最小值低至0，说明当波磨深度为240 μm时已经出现了轮轨分离现象，危及行车安全；整体来看不同波深对应轮轨力变化趋势一致，轮轨力与波深呈现正相关的关系.

图  5  不同波深轮轨力时程曲线

Figure  5.  Time history curves of wheel–rail force with different corrugation depths

下载: 全尺寸图片 幻灯片

值得注意的是，轮轨力时程曲线中给出了传统线弹性扣件模型波深0.12 m的轮轨竖向力时程曲线，对比非线性模型及线弹性模型轮轨竖向力变化清晰看出轮轨力对扣件频变特征的敏感性（如图5（a）中的红色圆圈所示）. 根据文献[1]中试验表明：扣件刚度与频率双对数函数成正比，当车轮行驶至扣件前端时受到轮轨振动频率变化影响，橡胶垫板刚度表现出刚硬特点使轮轨竖向力增大；当车轮运行至扣件正上方，在轴重作用下轮轨振动频率降低，橡胶垫板发生黏滞性流动，刚度减小，橡胶垫板表现出柔软特点使轮轨力减小，轮轨力出现先增大后减小的波动变化，2种扣件模型对轮轨力仿真结果影响偏差最大可达到13.1%.

图6为不同速度（250、300、350 km/h）工况下考虑橡胶垫板刚度非线性变化的仿真结果. 由图6可知：不同运行速度下轮轨力具有相似变化趋势和幅值，值得注意的是在0.5～1.0 m动态松弛区位置有明显的轮轨力滞后现象，即相同位置下速度越低轮轨力越大，而非线弹性模型仿真中相同位置轮轨力仅是幅值变化. 笔者认为发生此种现象原因为：轮轨动态响应触发非线性橡胶垫板刚度频变特性，不同的运行速度使得不同工况下车轮行驶至相同扣件位置所需时间不同，运行速度越小，给予非线性扣件橡胶垫板黏滞性流动的时间越长，橡胶垫板发生黏滞性流动导致扣件刚度增加，表现在轮轨力时程曲线上就是相同位置运行速度越小，轮轨力反而越大；换言之，运行速度越快，非线性扣件橡胶垫板来不及发生瞬态变化，车轮已完成行驶，相对而言速度快的工况轮轨力更小，造成了0.5～1.0 m动态松弛区不同运行速度轮轨力依次递减的现象；在波磨区并未体现出该种规律，这是因为波磨区域所激发的轮轨力达到了较大幅值，非线性扣件系统导致的轮轨力依次递减现象被削弱从而使得轮轨力时程曲线更加贴合而非依次递减.

图  6  不同运行速度轮轨力时程曲线

Figure  6.  Time history curves of wheel–rail force at different operating speeds

下载: 全尺寸图片 幻灯片

3.1.2   频域内轮轨力结果分析

本文基于扣件橡胶垫板刚度非线性变化特性建立车轮-轨道柔性有限元模型，通过对时域内提取数据进行小波变换得到轮轨力在频率影响下的计算结果是必要的. 对120、180、240 μm不同波磨波深以及250、300、350 km/h不同速度工况下频域结果进行对比分析，拾取轮轨力时程曲线中轮轨力变化显著的频段进行细致分析.

小波变换不同于Fourier变换，利用小波基伸缩和平移对时域内局部信号进行多尺度细化分析，根据小波支集长度、消失矩阶数以及正交性和对称性等小波选取原则^[14]，选择Gauss包络下的Morlet小波进行轮轨力时频域计算结果转换，因为Morlet小波具有较高的振幅聚集性^[15].

轮轨力时程曲线经过小波变换后得到的轮轨力频域响应如图7所示. 由图7可知：除初始激扰外高频响应主要出现在波深工况0.022～0.035 s位置即波磨求解区，并在0.030 s左右位置达到最大，频域响应图像呈规则棒状，图7（b） ② 显色卡深色分为2块是由于车轮进出波磨区产生2次明显频率变化，随波磨深度及速度增加，扣件胶垫受高频激扰而使得橡胶垫板刚度增加，进而加剧了轮轨冲击响应使显色卡着色加深.

图  7  不同波深/运行速度轮轨力时频图

Figure  7.  Time-frequency diagrams of wheel–rail force at different corrugation depths/operating speeds

下载: 全尺寸图片 幻灯片

图8为波磨区扣件位置处轮轨力的输出频谱图，对比分析不同速度和不同波深的主频响应及轮轨力变化. 从图8可知：随波深增加，轮轨激振响应的主频峰值位置相近，波磨深度越大对应频带越宽，分别对应858、882、904 Hz，对应轮轨力强度显著增加；第二峰值区的中心峰值处轮轨力同样存在小幅度增加；另外，不同运行速度对列车通过橡胶垫板时轮轨力幅值影响不大，但峰值主频出现明显变化，其对应频率分别为854、1025、1196 Hz.

图  8  不同波深/运行速度轮轨力功率谱曲线

Figure  8.  Power spectrum curves of wheel–rail force at different corrugation depths/operating speeds

下载: 全尺寸图片 幻灯片

3.2   橡胶垫板非线性力学特性对轴箱加速度影响

3.2.1   时域内轴箱加速度结果分析

波磨不平顺是影响列车行驶安全性、乘客乘车舒适性以及控制行车速度的主要原因之一，轴箱加速度是评价轨道几何不平顺，指导轮轨打磨旋修的重要指标^[16]，因此，相对车体及构架，基于轴箱加速度的诊断分析对轨道波磨伤损的预测更加科学有效.

图9为理想圆顺的车轮在无其他外部激扰的情况下不同波磨深度和运行速度工况轴箱加速度变化时程曲线. 由图9可知：由于施加速度和转矩，不可避免地于初始位置引入激扰导致轴箱加速度幅值变化，在动力松弛区有效衰减并保持平稳；进入波磨求解区后由于波磨的存在出现轴箱加速度波浪形规律变化，幅值变化与波深呈正相关，高速运行的列车通过波磨伤损时将引起更大的轮轨响应；当运行速度达到350 km/h时波磨区轴箱加速度变化幅值反而小于较低运行速度的轴箱加速度，笔者认为运行速度达到一定限值时轮轨产生不稳定接触，即车轮与波磨钢轨发生弹性接触使轮轨响应减弱，轴箱加速度相应降低. 对照3.1.1节轮轨力时域结果，350 km/h工况多次出现轮轨分离现象，高速列车驶过波磨区段存在严重的脱轨风险，应当按时检测车轮与钢轨服役状态并低速通过伤损区段，此外扣件位置车体轴箱加速度均显著降低，说明扣件系统橡胶垫板对降低列车振动、提高乘客乘车舒适性具有显著的作用.

图  9  不同波深/运行速度轴箱加速度时程曲线

Figure  9.  Time history curves of axle box acceleration at different corrugation depths/operating speeds

下载: 全尺寸图片 幻灯片

3.2.2   频域内轴箱加速度结果分析

图10为轴箱加速时频图，由图10可知：橡胶垫板非线性力学特性对于低频段（<200 Hz）没有明显的区别，而在高频段（800～1200 Hz）出现较大的轴箱加速度，显色卡的刻度反映出与试验测试结果相同的规律，轴箱加速度与波磨深度呈现正相关的关系；观察图10（b）显色卡刻度极值发现，运行速度提升引起波磨区轮轨振动更加剧烈，但当速度达到一定限值，由波磨引起的振动减弱，轴箱加速度降低.

图  10  轴箱加速度时频图

Figure  10.  Time-frequency diagrams of axle box acceleration

下载: 全尺寸图片 幻灯片

提取不同工况响应最剧烈时刻数据信息绘制功率谱曲线（图11），不同波深主频相近，频带随波深增大而变宽，相较而言频带宽度对运行速度更加敏感，250～300 km/h时主频增大与频带加宽显著.

图  11  不同波深/运行速度轴箱加速度功率谱曲线

Figure  11.  Power spectrum curves of axle box acceleration at different corrugation depths/operating speeds

下载: 全尺寸图片 幻灯片

4.   结　论

本文对高速铁路无砟轨道WJ-8扣件弹性橡胶垫板刚度进行非线性力学特性试验，提取非线性参数，运用ANSYS/Ls-dyna建立考虑橡胶垫板非线性力学特性的三维高速轮轨瞬态滚动接触有限元模型，再现扣件橡胶垫板真实服役状态刚度/阻尼非线性变化，探索基于非线性橡胶垫板的钢轨波磨处轮轨瞬态滚动接触动力相互作用，研究结论如下：

1） 轮轨力对扣件频变特性反应显著，相较线弹性扣件，在扣件位置处受到非线性扣件黏弹性特征变化影响，真实轮轨相互作用在扣件前端增大到达扣件上方后减小的波动形变化，2种扣件模型对轮轨力仿真结果最大差异可达到13.1%.

2） 不同运行速度导致轮轨力出现滞后现象，由于运行速度导致运行时间不同，非线性扣件橡胶垫板黏滞性流动的时间越长，扣件刚度增长越大造成0.5～1.0 m区段轮轨力与速度成反比现象.

3） 轮轨力时程曲线经Morlet小波变换绘制轮轨力频域响应，波磨求解区频域响应图像呈规则棒状，车轮进出波磨区产生两次明显频率变化引起扣件刚度变化致使速度工况时频图显色卡深色区域分为两块，分列波磨求解区两侧.

4） 提取车体轴箱加速度信号后发现，信号幅值和频带对运行速度变化更加敏感，250～300 km/h时主频增大与频带加宽显著，速度达到一定限值后信号减弱，300 km/h轴箱加速度信号最大达到119 m/s²，350 km/h最大主频为1044 Hz.

5） 本文所建立的扣件橡胶垫板刚度阻尼的非线性数据来源于试验，模型建立方法将橡胶垫板黏弹性频变特征纳入其中使得轮轨瞬态滚动接触分析更加准确可靠，推动了瞬态动力学研究中扣件橡胶垫板线弹性特征向非线性特征的转变，使得模型更加贴近实际，对今后瞬态动力学模型建立及仿真分析提供指导性建议和理论、试验依据.