基于SPH方法铁路车轴遭受道砟撞击的数值模拟*

敬 霖,韩亮亮,周彭滔

(西南交通大学牵引动力国家重点实验室,四川 成都 610031)

车轴是铁路机车、车辆走行部轮轴系统的重要组成部分，主要承受旋转弯曲载荷和扭转载荷作用。高速列车运行过程中，轨道上的道砟容易被卷起导致车轴轴身表面的冲击损伤，尤其是在雪季、隧道中和隧道重修后车轴可能遭受轨道之间的离散道砟和道床上的堆积道砟或其他尖锐硬物的撞击损伤[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个基本方程如下：

(1)

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

2 计算模型

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

高速轨道特级道砟实物如图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的正四面体道砟计算分析模型。

车轴材料为EA4T钢，采用塑性随动强化模型*MAT_PLASTIC_KINEMATIC描述，密度为7 850 kg/m3，杨氏模量为206 GPa，泊松比为0.3，切线模量为20.6 GPa，屈服应力为560 MPa[18]。道砟材质为花岗岩，是一种非均质多相脆性材料，这里采用*MAT_JOHNSON_HOLMQUIST_CERAMICS模型来描述，详细材料参数[19]：密度为2 657 kg/m3,剪切模量为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,D1=0.005,D2=0.7,K1=55.6 GPa,K2=-23 GPa,K3=2 980 GPa

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

(2)

通过定义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)波动直至整个撞击响应过程结束。

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种撞击速度下均近似等于零)。

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

3.2.2道砟尺寸的影响

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

3.2.3撞击角度的影响

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

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

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

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

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

Fmax=kδ+m,Fmax>Fcrit

(3)

式中:Fmax和δ分别为撞击力峰值和最大残余变形；Fcrit为产生残余变形的撞击力临界值(即只有当撞击力峰值大于Fcrit时车轴才会出现残余变形)，且有Fcrit=62.7 kN(正四面体道砟)和Fcrit=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可以看出，无论是正四面体道砟还是球形道砟撞击，给定撞击角度时，撞击力峰值与最大残余变形呈现出近似线性增大的关系，且随着撞击角度的增大，撞击力峰值随残余变形的增长率逐渐减小。也就是说，在相同撞击力载荷下，道砟与车轴的撞击角度越大，其残余变形就越大，撞击损伤就越严重。

3.4 车轴撞击损伤评估分析

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

参考文献：

[1] GRAVIER N, VIET J J, LELUAN A. Predicting the life of railway vehicle axles[C]∥Proceedings of the 12th International Wheelset Congress. Qingdao, China, 1998:133-146.

[2] 周素霞,谢基龙.高速客车空心轴车轴裂纹扩展特性研究[J].工程力学,2009,26(7):232-237.

ZHOU Suxia, XIE Jilong. Research of the fatigue crack propagation characteristic on railway hollow axles[J]. Engineering Mechanics, 2009,26(7):232-237.

[3] BERETTA S, GHIDINI A, LOMBARDO F. Fracture mechanisms and scale effects in the fatigue of railway axles[J]. Engineering Fracture Mechanics, 2005,72:195-208.

[4] OGNJANOVIC M, SIMONOVIC A, RISTIVOJEVIC M, et al. Research of rail traction shafts and axles fractures towards impact of service conditions and fatigue damage accumulation[J]. Engineering Failure Analysis, 2010,17:1560-1571.

[5] 刘志明,孙守光,缪龙秀.车轴裂纹扩展寿命的分析与计算方法[J].中国铁道科学,2008,29(3):89-94.

LIU Zhiming, SUN Shouguang, MIAO Longxiu. Analysis and calculation method of axle crack growth life[J]. China Railway Science, 2008,29(3):89-94.

[6] GOMEZ M J, CASTEJON C, GARCIA-PRADA J C. New stopping criteria for crack detection during fatigue tests of railway axles[J]. Engineering Failure Analysis, 2015,56:530-537.

[7] MCQUAID J, JONES N. A re-examination of Andrews’ research on impact resistance of railway axles[J]. International Journal of Impact Engineering, 1999,22:727-738.

[8] ANDREWS T. Effect of temperature on the strength of wrought iron railway axles: Part 2[J]. Minutes of the Proceedings, 1888,94:180-209.

[9] 周素霞.高速列车空心车轴损伤容限理论与方法研究[D].北京:北京交通大学,2009:71-80.

[10] 张俊清,周素霞,谢基龙.缺口对车轴钢疲劳性能的影响[J].北京交通大学学报,2010,34(4):133-135.

ZHANG Junqing, ZHOU Suxia, XIE Jilong. Experimental research on the effects of notches on the fatigue behavior of axle steel[J]. Journal of Beijing Jiaotong University, 2010,34(4):133-135.

[11] 刘军,李玉龙,刘元镛.基于SPH方法的叶片鸟撞数值模拟研究[J].振动与冲击,2008,27(9):91-93.

LIU Jun, LI Yulong, LIU Yuanyong. Numerical simulation study of bird-impact on a blade using SPH method[J]. Journal of Vibration and Shock, 2008,27(9):91-93.

[12] 李志强,韩强,杨建林,等.基于SPH方法的鸟撞飞机风挡的数值模拟[J].华南理工大学学报,2009,37(12):147-150.

LI Zhiqiang, HAN Qiang, YANG Jianlin, et al. SPH-based numerical simulation of aircraft windshield under bird impact[J]. Journal of South China University of Technology, 2009,37(12):147-150.

[13] 刘军,李玉龙,徐绯.基于PAM-CRASH的鸟撞飞机风挡动响应分析[J].爆炸与冲击,2009,29(1):80-84.

LIU Jun, LI Yulong, XU Fei. Dynamic response analysis of bird-impact aircraft windshields based on PAM-CRASH[J]. Explosion and Shock Waves, 2009,29(1):80-84.

[14] LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. Astronomical Journal, 1977,82(12):1013-1024.

[15] 孙晓旺,章杰,王肖钧,等.应力波数值计算中的SPH方法[J].爆炸与冲击,2017,37(1):10-14.

SUN Xiaowang, ZHANG Jie, WANG Xiaojun, et al. Application of SPH in stress wave simulation[J]. Explosion and Shock Waves, 2017,37(1):10-14.

[16] 铁道部运输局铁道科学研究院金属及化学研究所.铁路货车轮轴典型伤损图册[M].北京:中国铁道出版社,2006:64-66.

[17] 张徐,赵春发,翟婉明.铁路碎石道砟静态压碎行为数值模拟[J].西南交通大学学报,2015,50(1):137-143.

ZHANG Xu, ZHAO Chunfa, ZHAI Wanming. Numerical analysis of static crushed behavior of railway ballast[J]. Journal of Southwest Jiaotong University, 2015,50(1):137-143.

[18] 王筑生,梁益龙,吴少斌,等.大截面EA4T车轴钢回火工艺对组织和性能的影响[J].材料热处理学报,2012,33(5):48-52.

WANG Zhusheng, LIANG Yilong, WU Shaobin, et al. Effect of temperature process on microstructure and properties of EA4T axle steel[J]. Transactions of Materials and Heat Treatment, 2012,33(5):48-52.

[19] AI H A, AHRENS T J. Simulation of dynamic response of granite: A numerical approach of shock-induced damage beneath impact craters[J]. International Journal of Impact Engineering, 2006,33(1):1-10.

[20] LI Jun, WU Chengqing, HAO Hong. Investigation of ultra-high performance concrete slab and normal strength concrete slab under contact explosion[J]. Engineering Structures, 2015,102:395-408.

[21] HALLQUIST J O. LS-DYNA theory manual: Version 971[M]. Livermore: Livermore Software Technology Corporation, 2007:142-143.