基于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.