车轴是铁路机车、车辆走行部轮轴系统的重要组成部分，主要承受旋转弯曲载荷和扭转载荷作用。高速列车运行过程中，轨道上的道砟容易被卷起导致车轴轴身表面的冲击损伤，尤其是在雪季、隧道中和隧道重修后车轴可能遭受轨道之间的离散道砟和道床上的堆积道砟或其他尖锐硬物的撞击损伤^[1-2]。这些损伤不断积累和演化，可能导致车轴上疲劳裂纹的萌生和扩展，从而给列车的运行安全带来极大的隐患，严重时会造成车轴的疲劳断裂(冷切轴断裂)导致灾难性的重大倾覆事故^[3-4]。

目前，大多研究主要集中在铁路车轴的疲劳强度、裂纹检测及扩展寿命等方面^[2-6]，而针对车轴撞击损伤的研究较少。但是，铁路车轴抵抗撞击载荷的能力却是车轴早期设计中的一个主要考虑因素。由于当时还没有计及动态行为的车轴设计理论基础，一些铁路公司制定了自己的车轴撞击实验方法，最常用的方法是简支车轴跨中位置的落锤冲击实验^[7]。Andrews^[8]于1879年汇总了1870~1879年间的铸铁车轴撞击实验结果，实验中落锤的质量介于816.48~1 074.58 kg之间，下落高度位于5.334~9.144 m之间，冲击能量介于46.4~91.1 k J之间。McQuaid等^[7]基于Andrews整理的车轴撞击实验结果^[8]，将车轴简化为等截面的简支梁，给出了撞击载荷下车轴的最大残余变形的理论预测公式。然而，铁路车轴在服役过程中更容易遭受道砟等外物的撞击作用而产生撞击损伤，这些损伤与车轴材料的疲劳损伤、固有缺陷以及恶劣服役环境等因素相互耦合，势必会加剧裂纹的萌生和扩展。Gravier等^[1]较早提出了车轴受外物撞击擦伤的概念，他们指出这种撞击一旦发生，可能对车轴造成不同程度的损伤和尖锐缺口。周素霞^[9]将车轴受道砟冲击产生的缺口形状简化为三维半椭球形，开展了非动力车轴疲劳裂纹形成寿命预测。张俊清等^[10]采用有限元方法研究了车轴缺口形式对其材料疲劳寿命的影响。但是上述研究均未直接涉及车轴遭受道砟撞击的实验和有限元仿真分析。此外，道砟为脆性材料，在撞击车轴数值模拟过程中呈现大变形，常用的Lagrange算法难以避免有限元网格的严重畸变，影响数值计算精度，且对于高速撞击网格畸变会引起显式时间积分步长过短，极大提高了计算成本。采用ALE(arbitrary Lagrange Euler)算法虽然可以避免网格的严重畸变，但计算过程中需要不断地重构网格，从而导致计算时间的大幅度增加和计算精度的降低。SPH(smoothed particle hydrodynamics)算法既可避免Lagrange算法中大变形时网格畸变造成的精度破坏问题，还可避免Euler算法中Euler网格与材料的界面问题，在求解高速碰撞等动态大变形问题方面得到了广泛的应用^[11-13]。

本文中，基于SPH方法，采用非线性动力学有限元软件LS-DYNA模拟动车组车轴遭受道砟撞击时的动态力学响应，研究道砟几何形状、道砟尺寸、撞击速度和撞击角度对车轴撞击响应和损伤的影响。

1.   SPH算法

SPH方法是Lucy等^[14]于1977年提出的一种纯Lagrange、自适应性、无网格粒子法，最初主要用于解决三维开放空间的天体物理学问题。其中心思想是将问题域离散化为有限数量的粒子，每个粒子都具有独立的质量、密度及其他物理属性。随着计算方法的深入研究和实际应用的需求，SPH方法已被广泛应用于材料动态响应、具有大变形的流体动力学、侵彻和高速碰撞等方面。

SPH算法求解动力学问题主要包括2个方面：一是使用积分表达式对函数进行近似；二是通过对最近相邻粒子的值进行累加求和来近似离散点的函数值。其2个基本方程如下：

式中:f 为三维坐标变量x的任意函数； < f(x) >为函数f (x)的近似值；Ω为包含x的积分域；W(x－x′, h)被称为光滑函数或核函数，依赖于2点之间的距离|x－x′|和光滑长度h；i和j为粒子编号；m_j和ρ_j为粒子的质量和密度；N为粒子的总数；W_ij=W(x_i－x_j, h)。这里需要说明的是，不同的核函数形式对SPH算法仿真结果的影响较大，常见的核函数有B-样条函数、高斯型核函数和二次型核函数^[15]，本文中采用应用最广泛、也是LS-DYNA程序自嵌入的B-样条函数作为SPH算法的核函数。

2.   计算模型

以某动车组非动力车轴为研究对象，建立全尺寸车轴(长度为2 180 mm)的有限元模型，如图 1所示。采用八节点六面体实体单元，单元长度为5 mm，共487 424个单元。由于车轴受撞击擦伤的部位大多都集中在车轴中间部位^[16]，因此本文中选定车轴中间区域作为撞击响应的考察位置。

图  1  车轴的有限元模型

Figure  1.  The finite element model of the axle

下载: 全尺寸图片 幻灯片

高速轨道特级道砟实物如图 2(a)所示，直径范围一般为20~60 mm^[17]。本文中分别将道砟模型简化为球形和正四面体，采用等效体积法，选定20、40和60 mm等3种直径规格和33、66和99 mm等3种边长尺寸，建立了不同规格道砟的SPH模型。3种直径球形道砟SPH模型的粒子数目分别为1 791、8 217和22 575，3种边长尺寸正四面体道砟模型的粒子数目分别为1 620、8 181和21 826。图 2(b)和(c)分别给出了直径为40 mm的球形道砟和边长为66 mm的正四面体道砟计算分析模型。

图  2  道砟实物和SPH模型

Figure  2.  Photo and computational models of the ballast

下载: 全尺寸图片 幻灯片

车轴材料为EA4T钢，采用塑性随动强化模型*MAT_PLASTIC_KINEMATIC描述，密度为7 850 kg/m³，杨氏模量为206 GPa，泊松比为0.3，切线模量为20.6 GPa，屈服应力为560 MPa^[18]。道砟材质为花岗岩，是一种非均质多相脆性材料，这里采用*MAT_JOHNSON_HOLMQUIST_CERAMICS模型来描述，详细材料参数^[19]：密度为2 657 kg/m³, 剪切模量为30 GPa，参考应变率为0.001 s^-1，拉伸强度为0.15 GPa, 最大断裂强度为0.2 GPa，Hugoniot弹性极限强度为2.66 GPa，Hugoniot弹性极限压力为2.73 GPa；其他参数a=1.01，b=0.68，c=0.005，M=0.76，N=0.83，β=1.0, D₁=0.005, D₂=0.7, K₁=55.6 GPa, K₂=-23 GPa, K₃=2 980 GPa

道砟状态方程选用*EOS_GRÜNEISEN，具体参数^[20]为：C₀=2 100 m/s, S₁=1.4, γ₀=2。*MAT_JOHNSON_HOLMQUIST_CERAMICS模型中材料的量纲一等效应力^[21]表达式为：

式中：σ_i^*和σ_f^*分别为量纲一的等效应力和断裂强度；D为损伤系数；a和N为完整强度常数；p^*为量纲一压力；T^*为量纲一拉伸强度；${\dot \varepsilon ^*}$为量纲一塑性应变率；c为应变率系数；b和M为断裂强度常数。

通过定义SPH单元为从动面，车轴有限元模型为主动面，选用*CONTACT_ERODING_NODES_TO_SURFACE接触算法来实现有限元单元和SPH单元的耦合。对车轴两端面节点实施完全固支边界条件。

3.   模拟结果与分析

3.1   车轴遭受道砟撞击的响应过程

图 3给出了边长为99 mm的正四面体道砟以v=400 km/h的速度撞击车轴时的等效应力云图。为便于观察，图 3同时给出了车轴受道砟撞击区域的局部等效应力云图及其对应时刻车轴撞击截面的等效应力云图。从图中可以看出，道砟以给定速度撞击车轴时，车轴受到的等效应力在碰撞瞬间急剧增大，t =0.15 ms时车轴撞击位置出现最大的等效应力, 约583.5 MPa；随后向车轴内部和两侧发展，其等效应力逐渐减小至249.6 MPa(t =0.20 ms)；此后随着道砟SPH粒子与车轴的进一步耦合作用，车轴的等效应力稍有增大(约282.1 MPa)；最后在一定幅值范围内(275.4~277.8 MPa)波动直至整个撞击响应过程结束。

3a  t=0 ms时道砟撞击车轴的典型等效应力云图

3a.  Typical effective stress contours of the axle subjected to ballast impact at t=0 ms

下载: 全尺寸图片 幻灯片

3b  t=0.15 ms时道砟撞击车轴的典型等效应力云图

3b.  Typical effective stress contours of the axle subjected to ballast impact at t=0.15 ms

下载: 全尺寸图片 幻灯片

3c  t =0.20 ms时道砟撞击车轴的典型等效应力云图

3c.  Typical effective stress contours of the axle subjected to ballast impact t =0.20 ms

下载: 全尺寸图片 幻灯片

3d  t =0.27 ms时道砟撞击车轴的典型等效应力云图

3d.  Typical effective stress contours of the axle subjected to ballast impact at t =0.27 ms

下载: 全尺寸图片 幻灯片

3e  t =0.40 ms时道砟撞击车轴的典型等效应力云图

3e.  Typical effective stress contours of the axle subjected to ballast impact at t =0.40 ms

下载: 全尺寸图片 幻灯片

3f  t =0.60 ms时道砟撞击车轴的典型等效应力云图

3f.  Typical effective stress contours of the axle subjected to ballast impact at t =0.60 ms

下载: 全尺寸图片 幻灯片

3.2   撞击力和变形响应特征

撞击力和变形是道砟撞击车轴分析与评价中最关键的2个指标。本文中比较了2种形状道砟在不同撞击速度、道砟尺寸和撞击角度工况下的撞击力响应和变形响应。

3.2.1   撞击速度的影响

考虑到动车组的运行速度及列车交会等实际服役工况，仿真计算中分别选取v=100，200，300，400 km/h等4种工况进行讨论。图 4给出了3种边长的正四面体道砟以不同速度撞击车轴时的力-时间响应曲线以及撞击点的位移-时间响应曲线。可以看出，同一边长条件下，随着撞击速度的增大，撞击力峰值明显增大，而撞击响应时间先逐渐增大后趋于稳定(当边长为33和66 mm时，撞击响应时间在撞击速度大于等于300 km/h时趋于稳定，而当边长为99 mm时，撞击响应时间在撞击速度大于等于200 km/h时趋于稳定)。并且，随着撞击速度的增大，撞击力-时程曲线呈现出明显的“锯齿形”振荡。这是由于撞击速度较高时，道砟的初始动量很难在短时间内衰减至零，道砟SPH粒子与车轴的多次耦合作用造成的。从图 4还可以看出，车轴受撞击区域中心点的最大瞬态变形和残余变形均随着撞击速度的增大而增大(道砟边长为33 mm时，车轴受撞击点的残余变形在4种撞击速度下均近似等于零)。

图  4  不同速度下正四面体道砟撞击车轴时的力与变形响应时程曲线

Figure  4.  Force and displacement history curves of the axle subjected to the regular-tetrahedron-ballast impact at different velocities

下载: 全尺寸图片 幻灯片

为了比较，图 5给出了3种直径球形道砟在不同撞击速度下车轴的力-时间和位移-时间响应曲线。同样地，随着撞击速度的增大，撞击力峰值增大，撞击响应时间也相应增大，且响应时间在撞击速度大于等于300 km/h时也趋于一致；车轴受撞击区域中心点的最大瞬态变形和残余变形(除直径为20 mm时残余变形近似等于零之外)均随着撞击速度的增大而相应增大。结合图 4和图 5可以看出，3种直径球形道砟在不同速度撞击下的撞击力峰值均大于对应的四面体道砟撞击产生的撞击力峰值，且振荡明显减小；车轴受撞击点变形在道砟尺寸较小时(边长小于等于66 mm或直径小于等于40 mm)在不同速度下受道砟形状的影响较小，但当道砟边长为99 mm(或直径60 mm)时，不同速度下四面体道砟撞击车轴产生的最大瞬态变形和残余变形均大于相应的球形道砟工况，这是因为道砟尺寸较大时，道砟的初始动能也就较大，四面体道砟以尖锐棱边撞击车轴时产生了较大的局部塑性变形。

图  5  不同速度下球形道砟撞击车轴时的力与变形响应时程曲线

Figure  5.  Force and displacement history curves of the axle subjected to the spherical ballast impact at different velocities

下载: 全尺寸图片 幻灯片

3.2.2   道砟尺寸的影响

铁路车轴服役过程中，被卷起的道砟尺寸不同时对车轴的撞击损伤也必然不同，因此有必要分析不同尺寸的道砟对车轴的撞击响应。明显地，图 4和图 5中已经包括了3种尺寸规格道砟对车轴撞击响应的影响。为了简洁，这里仅分析速度为400 km/h时道砟尺寸的影响，如图 6所示。给定撞击速度时，道砟尺寸越大，其初始动能就越大，撞击力峰值及响应时间也相应增大，撞击车轴后产生的瞬态和残余变形也就越大。如上所述，球形道砟撞击车轴产生的撞击力峰值均大于对应等效尺寸的正四面体道砟撞击产生的撞击力峰值，但残余变形却小于正四面体道砟撞击的情形。

图  6  不同尺寸道砟以400 km/h的速度撞击车轴时的力与变形响应时程曲线

Figure  6.  Force and displacement history curves of the axle subjected to the impact of the ballastswith different sizes at 400 km/h

下载: 全尺寸图片 幻灯片

3.2.3   撞击角度的影响

实际情况中，道砟与车轴的碰撞行为并不一定是理想的垂直正碰，可能会出现道砟与车轴中心线成一定角度的斜碰撞。本文中比较了道砟与车轴中心线分别成45°、60°和90°等3种碰撞工况下车轴的撞击响应。图 7给出了正四面体道砟以不同撞击速度、不同角度撞击车轴时的力和位移响应时程曲线。从图中可以看出，对于4种撞击速度工况，不同撞击角度下车轴的力和位移响应曲线变化趋势大致相同，且随着撞击角度的增大，撞击力峰值相应增大，车轴撞击点的最大瞬态变形和最终残余变形也随之增大。显然，3种撞击角度工况中，道砟与车轴中心线成90°垂直正碰时车轴受到的撞击力和变形均为最大，属于碰撞行为中最危险的工况。

7a  正四面体道砟以100 km/h的速度、不同的角度撞击车轴时的力和位移响应曲线

7a.  Force and displacement response history curves of the axle subjected to the regular-tetrahedron ballast impact with different angles at 100 km/h

下载: 全尺寸图片 幻灯片

7b  正四面体道砟以200 km/h的速度、不同的角度撞击车轴时的力和位移响应曲线

7b.  Force and displacement response history curves of the axle subjected to the regular-tetrahedron ballast impact with different angles at 200 km/h

下载: 全尺寸图片 幻灯片

7c  正四面体道砟以300 km/h的速度、不同的角度撞击车轴时的力和位移响应曲线

7c.  Force and displacement response history curves of the axle subjected to the regular-tetrahedron ballast impact with different angles at 300 km/h

下载: 全尺寸图片 幻灯片

7d  正四面体道砟以400 km/h的速度、不同的角度撞击车轴时的力和位移响应曲线

7d.  Force and displacement response history curves of the axle subjected to the regular-tetrahedron ballast impact with different angles at 400 km/h

下载: 全尺寸图片 幻灯片

图 8给出了球形道砟以不同初始速度和角度撞击车轴时的撞击力和位移响应时程曲线。与正四面体道砟撞击一样，无论在哪种速度工况下，撞击力峰值、车轴最大瞬态变形和残余变形均随着撞击角度的增大而增大，道砟与车轴中心线成90°的垂直正碰是最危险的工况。相同撞击速度下，球形道砟以3种不同角度(45°、60°和90°)撞击车轴产生的撞击力峰值均大于对应的正四面体道砟撞击情况，但其车轴的残余变形均小于对应的正四面体道砟撞击情况，这是由不同的变形模式和作用机理造成的。

图  8  球形道砟以不同初始速度和角度撞击车轴时的撞击力和位移响应时程曲线

Figure  8.  Force and displacement response history curves of the axle subjected to the spherical ballast impact with different angles at different impact velocities

下载: 全尺寸图片 幻灯片

3.3   车轴残余变形与峰值力之间的关系

图 9给出了正四面体和球形道砟与车轴中心线成90°垂直正碰车轴时撞击力峰值与最大残余变形之间的关系。可以看出，当撞击载荷峰值低于一定值时，车轴变形处于弹性变形范围内，一旦撞击力峰值超过其极限值，车轴将会出现残余变形，且最大残余变形值与撞击力峰值呈近似线性增大的关系。球形道砟撞击时，撞击力峰值迅速增大，但相应的残余变形增量并不明显; 而正四面体道砟撞击时，随着撞击力峰值的缓慢增大，相应的残余变形则显著增大。也就是说，相同条件下车轴抵抗球形道砟撞击的性能优越于抵抗正四面体道砟撞击的性能。这是因为道砟撞击车轴时，车轴吸收的冲击能以塑性变形能的形式耗散，正四面体道砟撞击时车轴轴身的局部塑性变形较大，且撞击响应时间也比球形道砟撞击时长。

图  9  车轴遭受道砟撞击残余变形与撞击力峰值的关系

Figure  9.  Relationship between peak impact force and residual deformation of the axle subjected to ballast impact

下载: 全尺寸图片 幻灯片

通过数据拟合，可得到车轴遭受道砟撞击时最大残余变形与撞击力峰值之间的定量关系：

式中:F_max和δ分别为撞击力峰值和最大残余变形；F_crit为产生残余变形的撞击力临界值(即只有当撞击力峰值大于F_crit时车轴才会出现残余变形)，且有F_crit=62.7 kN(正四面体道砟)和F_crit=126.6 kN(球形道砟)；k和m为拟合参数，且对正四面体道砟有k=779.2 kN/mm和m=64.6 kN，对球形道砟有k=5 255 kN/mm和m=122.2 kN。

图 10(a)和(b)分别给出了正四面体和球形道砟以不同角度撞击车轴时撞击力峰值与最大残余变形之间的关系。从图 10可以看出，无论是正四面体道砟还是球形道砟撞击，给定撞击角度时，撞击力峰值与最大残余变形呈现出近似线性增大的关系，且随着撞击角度的增大，撞击力峰值随残余变形的增长率逐渐减小。也就是说，在相同撞击力载荷下，道砟与车轴的撞击角度越大，其残余变形就越大，撞击损伤就越严重。

图  10  不同撞击角度下车轴遭受道砟撞击残余变形与撞击力峰值的关系

Figure  10.  Relationship between peak impact force and residual deformation of the axle subjected to ballast impact at different angles

下载: 全尺寸图片 幻灯片

3.4   车轴撞击损伤评估分析

这里，采用车轴撞击部位的最大残余变形(压痕深度)来描述其撞击损伤，分析车轴最大压痕深度δ与吸收冲击能E的关系。图 11给出了道砟垂直正碰撞车轴时量纲一压痕深度δ与吸收冲击能平方根E^1/2之间的关系。从图中可以看出，车轴吸收的冲击能较小时，残余变形几乎为零，不会产生撞击损伤，而当吸收冲击能大于某一临界值时，量纲一压痕深度δ与吸收冲击能平方根E^1/2呈近似线性关系。通过数据拟合，可得到两者的函数关系式：

图  11  车轴量纲一压痕深度与吸收冲击能平方根之间的关系

Figure  11.  Relationship between the normalized indentation depth and the square root of impact energy

下载: 全尺寸图片 幻灯片

式中:α和β为拟合参数。对于正四面体道砟, 有δ=δ/L，E_crit = 0.24 J，α=5.47×10^-4 J^-1/2和β=-2.66×10^-4；对于球形道砟, 有δ=δ/D，E_crit = 0.15 J，α=4.28×10^-4 J^-1/2和β =-1.49×10^-4。

从图 4~8可以得到，铁路车轴在服役过程中遭受道砟撞击时轴身可产生最大残余变形(压痕深度)为0.12 mm的撞击坑。据相关统计，高速列车在有砟轨道运行时车轴遭受道砟撞击轴身表面可产生深度为0.3~0.8 mm的撞击坑。本文仿真结果略小于实测统计数据，这是因为仿真分析中车轴为理想的未服役车轴，一定程度上强化了车轴的抵抗道砟撞击的性能。如果考虑车轴材料的固有缺陷(微裂纹、夹杂等)、疲劳损伤积累以及低温服役环境(如哈大线最低环境温度可达-40 ℃)等因素的耦合作用，相同条件下车轴的撞击损伤必定更严重。这些损伤在周期性轮轨接触载荷下，尤其当轮轨界面存在接触不平顺(如车轮不圆顺、轨道波磨和轨接头等)引起的轮轨冲击载荷作用下，势必会加剧车轴的疲劳裂纹萌生和扩展，从而给铁道车辆带来了安全隐患。

4.   结论

轨道道砟撞击产生的轴身表面损伤是诱使铁路车轴疲劳裂纹萌生和扩展的原因之一，道砟撞击车轴属于典型的接触-碰撞弹塑性变形动力学计算分析问题，该碰撞过程包含几何非线性、物理非线性与接触非线性等问题。本文中基于SPH方法的数值模拟研究，得到了如下主要结论：

(1) 道砟撞击铁路车轴是发生在毫秒量级的冲击力学行为，撞击响应时间约在0.5 ms以内。

(2) 撞击力峰值和车轴残余变形随着撞击速度、道砟尺寸和撞击角度的增大也相应增大，车轴最大残余变形值与撞击力峰值呈近似线性增大的关系。

(3) 道砟与车轴中心线成90°的垂直正碰是碰撞行为中最危险的工况，且正四面体道砟对车轴轴身的撞击损伤较球形道砟撞击更严重，可产生最大残余变形(撞击坑深度)约为0.12 mm。

(4) 道砟与车轴垂直正碰情况下，车轴的量纲一压痕深度与吸收冲击能平方根成线性关系。