于榮泉, 任尊松, 孫守光

(北京交通大學 軌道車輛結構可靠性與運用檢測技術教育部工程研究中心， 北京 100044)

輪軌滾動接觸疲勞破壞現(xiàn)象隨著鐵路客貨運量的增大和列車速度的提高變得越來越嚴重，目前已成為影響和制約輪軌服役壽命和運輸安全的重要因素，是輪軌的主要失效形式之一[1]。1998年德國漢堡開往慕尼黑的ICE高速列車因低噪聲彈性車輪崩裂而脫軌撞橋，造成嚴重脫軌事故。2000年，一列高速列車從倫敦駛往Leeds的途中，由于曲線外側(cè)鋼軌疲勞裂紋擴展導致鋼軌斷裂造成出軌事故。由此可見研究輪軌滾動接觸疲勞損傷問題以降低由滾動接觸疲勞造成的列車事故具有重要的工程意義。滾動接觸彈塑性應力應變分析是研究輪軌滾動接觸疲勞失效的基礎和關鍵。

目前國內(nèi)外學者利用有限元方法對滾動接觸彈塑性應力做了很多研究。文獻[2-3]利用有限元軟件ABAQUS建立非穩(wěn)態(tài)載荷作用下鋼軌彈塑性二維有限元模型，并分析了波狀表面波谷和波峰處的殘余應力和殘余應變；文獻[4]建立了輪軌接觸熱彈塑性平面應變熱機耦合二維有限元模型，分析摩擦生熱對輪軌熱彈塑性的影響；文獻[5]用有限元方法建立了鋼軌三維彈塑性滾動接觸計算模型，模型中采用雙線性隨動硬化彈塑性本構模型，分析了鋼軌材料屈服強度對鋼軌殘余應力和應變的影響；文獻[6]利用有限元參數(shù)二次規(guī)劃法進行不同輪徑、軸重和摩擦系數(shù)等工況下彈性和彈塑性計算，得出輪軌間接觸狀態(tài)和接觸內(nèi)力的分布情況；文獻[7]建立基于ALE有限元方法的穩(wěn)態(tài)輪軌滾動接觸的三維有限元模型，該模型在虛功率方程中通過Lagrange乘子法引入接觸界面上無切向滑移約束，計算分析接觸斑的黏著特性；文獻[8]利用有限元方法，采用非線性各向同性隨動硬化模型，對循環(huán)載荷作用下鋼軌材料的棘輪效應進行分析。以上研究僅局限于輪軌接觸中鋼軌的彈塑性應力應變及接觸狀態(tài)和黏著特性的分析，有關車輪三維彈塑性應力應變分析的文獻較少。文獻[9]研究了車輪三維滾動接觸狀態(tài)下等效塑性應變和靜水壓力等隨軸重和溫度的變化規(guī)律，但沒有分析車輪表面附近殘余應力和殘余應變分布特性及隨循環(huán)載荷的累積規(guī)律。

鑒于此，本文運用有限元方法建立車輪三維滾動接觸彈塑性有限元計算模型，模型中采用Lemaitre-Chaboche非線性各向同性隨動硬化循環(huán)塑性模型，研究在循環(huán)載荷作用下，最大接觸壓力和切向力對車輪接觸表面附近材料的殘余應力和殘余應變分布特性及累積規(guī)律的影響。

1 數(shù)值計算模型

1.1 三維滾動接觸

利用ABAQUS軟件建立車輪三維彈塑性有限元模型。模型中不考慮循環(huán)載荷作用下車輪接觸表面材料塑性變形對接觸邊界條件的影響，輪軌表面接觸壓力和接觸斑形狀大小依據(jù)Hertz接觸理論計算[10]。

Hertz接觸理論中接觸斑法向接觸壓力分布(以接觸斑中心為原點，車輪周向為x軸方向，軸向為y軸方向，見圖1。接觸壓力為

( 1 )

式中：px,y為橢圓接觸斑內(nèi)任意點的接觸壓力；x和y分別為接觸斑周向和軸向坐標；a和b為橢圓形接觸斑長短半軸，a、b具體值由輪軌接觸幾何參數(shù)和輪軌力確定[11]。

( 2 )

式中：p0為接觸斑最大接觸壓力；Fn為輪軌接觸力。

車輪橢圓接觸斑長半軸和短半軸大小為

( 3 )

( 4 )

式中：mH、nH為赫茲接觸系數(shù)，與(B-A)/(B+A)的值有關；A、B為相對曲率，具體表達式見文獻[11]；kH1和kH2是與輪軌材料相關的常數(shù)，為

( 5 )

( 6 )

式中：E1和E2分別為車輪和鋼軌材料的彈性模量，取為210 GPa；v1和v2分別為車輪和鋼軌材料的泊松比，取為0.3。

輪軌幾何尺寸分別為：車輪滾動圓半徑R11為420 mm，車輪踏面橫向圓弧半徑為無窮大，鋼軌縱向曲率半徑為無窮大，鋼軌橫向曲率半徑R22為300 mm。

輪軌接觸斑內(nèi)切向力分布不考慮接觸斑黏著區(qū)和滑動區(qū)的影響，假設切向力正比于法向壓力。

Qx=μP

( 7 )

式中：Qx為x向切向力；μ為全滑動摩擦因數(shù)；P為法向接觸壓力。

1.2 有限元網(wǎng)格及其邊界條件

通過非線性軟件ABAQUS建立車輪三維彈塑性有限元模型，為了提高分析效率和節(jié)約計算成本，選取車輪圓心角為15°的扇形區(qū)域作為車輪計算模型，見圖2(a)。通過試算，循環(huán)載荷對車輪材料彈塑性影響主要分布在車輪表面及表面下15 mm深度范圍內(nèi)，在該范圍內(nèi)進行網(wǎng)格細化處理，在車輪表面徑向深度15 mm，軸向?qū)挾?5 mm的范圍內(nèi)采用單元尺寸為1 mm的正六面體縮減積分單元，遠離接觸面的區(qū)域采用較粗的網(wǎng)格進行處理，網(wǎng)格劃分細節(jié)見圖2(b)。車輪扇形有限元計算模型，徑向側(cè)面邊界節(jié)點采用周向和軸向位移約束邊界，輪轂與車軸配合表面節(jié)點采用位移全約束邊界。

1.3 材料本構模型

材料本構模型采用Lemaitre-Chaboche非線性各向同性隨動硬化循環(huán)塑性模型[12]，該模型能夠合理地描述滾動接觸條件下材料循環(huán)塑性的應力松弛和棘輪效應。材料參數(shù)采用鐵路工業(yè)常用的BS11普通車輪鋼，車輪鋼材料參數(shù)參見表1。

表1 車輪鋼循環(huán)塑性模型材料參數(shù)[13-14]

1.4 模型計算參數(shù)

為了計算和表達方便，文中對長度、應力和應變變量分別用a0、k和k/G進行無量綱化處理，a0為接觸斑長軸之半，k和G分別為剪切屈服強度和剪切彈性模量。

車輪三維彈塑性有限元模型主要分析循環(huán)接觸載荷反復作用下最大接觸壓力p0和切向力Qx對車輪接觸表面及附近材料殘余應力和殘余應變分布特性及累積規(guī)律的影響。選取軸重21、25、30 t，摩擦因數(shù)0、0.1、0.3、0.5分別進行計算。

軸重21、25、30 t分別用最大接觸壓力p0=3.97k、4.21k、4.47k表示；摩擦因數(shù)0、0.1、0.3、0.5分別表示Qx為0、0.1P、0.3P、0.5P。由式( 3 )和式( 4 )得到，p0=4.21k時接觸斑長軸之半為a0=7.1 mm。對車輪接觸表面及附近材料殘余應力和殘余應變分布特性及累積規(guī)律的影響。選取軸重21、25、30 t，摩擦因數(shù)0、0.1、0.3、0.5分別進行計算。

2 結果分析

考慮到輪軌滾動接觸三維彈塑性有限元分析的困難性和計算量巨大，文中對每種載荷工況循環(huán)計算10次。

2.1 滾動接觸彈塑性分析

圖3和圖4分別給出了p0=4.21k、Qx=0.3P工況下，接觸斑中心深度為z=0.067a0處應力、應變在第一次循環(huán)載荷作用下隨載荷步的變化。從圖3、圖4中可以看出，經(jīng)歷第一次循環(huán)載荷作用后，車輪三維滾動接觸狀態(tài)下6個殘余應力分量和6個殘余應變分量均不為零，這與文獻[15]分析結果一致。本文重點分析周向殘余應力σxr、軸向殘余應力σyr、殘余剪切應變γxzr隨循環(huán)載荷的變化情況。

圖5給出了p0=4.21k、Qx=0.3P工況下，σxr、σyr隨深度分布的累積情況。由圖5中可以看出，最大殘余應力σxr和σyr分別發(fā)生在車輪次表面z/a0=0.6和z/a0=0.4處。周向和軸向殘余應力在接觸表面和次表面均為殘余壓應力。

圖5(a)中，當z/a0>1.1時，產(chǎn)生周向殘余拉應力；圖5(b)中，當z/a0>0.9時，產(chǎn)生軸向殘余拉應力。從圖5中可以看出，經(jīng)過10次循環(huán)載荷作用后，殘余應力逐漸趨于穩(wěn)定。

圖6給出了p0=4.21k、Qx=0.3P工況下，γxzr隨深度的分布累積情況。從圖中可以看出，最大殘余剪切應變發(fā)生在接觸表面，且殘余剪切應變主要分布在接觸表面至z/a0=1范圍內(nèi)。殘余剪切應變隨循環(huán)載荷作用次數(shù)的增加而逐漸增大，但增長率隨循環(huán)載荷作用次數(shù)的增加而逐漸衰減。

2.2 最大接觸壓力影響

最大接觸壓力分別為p0=3.97k、4.21k和4.47k時，循環(huán)載荷作用10次后σxr、σyr隨深度的分布情況見圖7。

圖7(a)和圖7(b)分別為σxr和σyr的分布情況。從圖中可以看出，最大接觸壓力對殘余應力分布的影響較大。周向和軸向最大殘余應力均發(fā)生在次表面z/a0=0.35處，且隨著最大接觸壓力的增加而增大。

圖8所示為循環(huán)載荷作用10次后最大接觸壓力對殘余剪切應變的影響。由此可見，最大接觸壓力作用下最大殘余剪切應變均發(fā)生在z/a0=0.5處，且隨著最大接觸壓力的增加而增大。

2.3 切向力影響

圖9(a)和圖9(b)分別給出了循環(huán)載荷p0=4.21k作用10次后，切向力對周向殘余應力σxr和軸向殘余應力σyr的影響。從圖9(a)中可以看出，切向力對周向殘余應力σxr分布的影響顯著。切向力為0P、0.1P和0.3P時，殘余應力σxr最大值均發(fā)生在次表面，且切向力為0P和0.1P時發(fā)生在次表面z/a0=0.36處，切向力為0.3P時發(fā)生在次表面z/a0=0.64處。當切向力為0.5P時，殘余應力σxr最大值發(fā)生在車輪接觸表面。

由圖9(b)可得到，相比于圖9(a)，在不同切向力作用下，軸向殘余應力σyr最大值均發(fā)生在次表面z/a0=0.36處。純滾動和Qx=0.1P時，接觸表面產(chǎn)生軸向殘余拉應力，且2種工況下軸向殘余應力隨深度的分布曲線基本重合。Qx=0.3P和Qx=0.5P時，接觸表面產(chǎn)生軸向殘余壓應力。

圖10所示為循環(huán)載荷p0=4.21k作用10次后，切向力對殘余剪切應變γxzr隨深度分布的影響?？梢钥闯?，切向力Qx對殘余剪切應變的影響較大。在接觸表面附近，γxzr隨切向力的增加以指數(shù)形式增長，與文獻[15]分析結果趨勢相似。當Qx/P大于等于0.3時，最大γxzr發(fā)生在接觸表面，與文獻[16]試驗觀察一致。當Qx/P=0時，次表面γxzr為負值，但相比其他3種工況，剪切應變很小。3種Qx/P非零工況，在當前坐標系下，殘余剪切應變?yōu)樨撝怠Ｔ?種Qx/P工況下，殘余剪切應變均分布在z/a<1.1范圍內(nèi)，Qx/P對殘余剪切應變深度分布范圍影響較小。

3 結論

利用ABAQUS建立了車輪三維彈塑性有限元模型，模型中采用Lemaitre-Chaboche非線性各向同性隨動硬化循環(huán)塑性模型，運用Hertz接觸理論和全滑動切向理論來模擬輪軌反復滾動過程，通過有限元分析得到以下結論：

(1) 通過數(shù)值模擬得到了車輪接觸表面附近的殘余應力和殘余應變隨循環(huán)載荷的累積發(fā)展過程和隨深度的分布規(guī)律。

(2) 純滾動工況下最大接觸壓力越大，殘余應力和殘余應變越大，殘余應力和殘余應變最大值均發(fā)生在車輪次表面。

(3) 切向力對車輪接觸表面附近的殘余應力和殘余應變影響顯著，切向力作用下殘余剪切應變最大值發(fā)生在滾動接觸表面。

參考文獻：

[1] 金學松，張繼業(yè)，溫澤峰，等. 輪軌滾動接觸疲勞現(xiàn)象分析[J]. 機械強度, 2002, 24(2):250-257.

JIN Xuesong，ZHANG Jiye，WEN Zefeng，et al. Overview of Phenomena of chanical Strength，2002，24(2):250-257.

[2] 溫澤峰，金學松. 非穩(wěn)態(tài)純滾動接觸彈塑性分析[J]. 固體力學學報，2006，27(4)：355-361.

WEN Zefeng，JIN Xuesong. Elastic-plastic Analysis of Non-steady State Pure Rolling Contact[J].Acta Mechanica Solida Sinica，2006，27(4)：355-361.

[3] 溫澤峰，金學松，肖新標.非穩(wěn)態(tài)載荷對二維輪軌純滾動接觸應力和變形的影響[J].交通運輸工程學報，2006，27(4)：14-19.

WEN Zefeng, JIN Xuesong, XIAO Xinbiao. Influence of Non-steady State Loading on Two-dimensional Wheel-rail Pure Rolling Contact Stresses and Deformation[J]. Journal of Traffic and Transportation Engineering, 2006,27(4)：14-19.

[4] 李偉，溫澤峰，吳磊，等. 車輪滑動時鋼軌熱彈塑性有限元分析[J]. 機械工程學報，2010，46(10)：95-101.

LI Wei, WEN Zefeng, WU Lei, et al.Thermo-elasto-plastic Finite Element Analysis of Rail During Wheel Sliding[J]. Journal of Mechanical Engineering, 2010, 46(10); 95-101.

[5] 郭俊，溫澤峰，金學松，等.鋼軌三維彈塑性滾動接觸應力[J]. 西南交通大學學報，2007，42(3)：262-268.

GUO Jun，WEN Zefeng，JIN Xuesong，et al. Three Dimensional Elastic-plastic Rolling Contact Stresses in Rail[J]. Journal of Southwest Jiao Tong University, 2007, 42(3)：262-268.

[6] 張軍，吳昌華.輪軌接觸問題的彈塑性分析[J].鐵道學報,2000, 22(3): 16-21.

ZHANG Jun, WU Changhua. Elasto-plastic Analysis of Wheel Rail Contact Problem[J]. Journal of the China Railway Society,2000,22(3): 16-21.

[7] 常崇義，王成國.基于ALE有限元的輪軌穩(wěn)態(tài)滾動接觸分析[J]. 中國鐵道科學，2009,30(2):87-93.

CHANG Chongyi, WANG Chengguo.Wheel-rail Steady State Rolling Contact Analysis Based on ALE Finite Element Method[J]. China Railway Science, 2009,30(2):87-93.

[8] 肖乾，車宇翔，周新建，等.輪軌滾動接觸棘輪效應數(shù)值分析[J].鐵道學報，2013,35(12):19-23.

XIAO Qian, CHE Yuxiang, ZHOU Xinjian, et al. Numerical Analysis on Ratcheting Effect of Rolling Contact Between Wheel and Rail[J].Journal of the China Railway Society,2013,35(12):19-23.

[9] VERNERSSON T,CAPRIOLI S,KABO E, et al. Wheel Tread Damage: A Numerical Study of Railway Wheel Tread Plasticity Under Thermomechanical Loading[J]. Part F: J. Rail and Rapid Transit，2010,244(5)：435-443.

[10] JOHNSON K L. Contact Mechanics[M]. Cambridge:Cambridge University Press，1985.

[11] 金學松，劉啟躍.輪軌摩擦學[M].北京：中國鐵道出版社，2004：51-57.

[12] LEMAITRE J, CHABOCHE J L. Mechanics of Solid Materials[M]. Cambridge:Cambridge University Press, 1990.

[13] EKH M, JOHANSSON A, THORBERNTSSON H, et al. Models For Cycilc Ratchetting Plasticity-Integration and Calibration[J].Journal of Engineering Materials and Technology,2000,122(6):49-55.

[14] RINGSBERG J W. Cyclic Ratchetting and Failure of a Pearlitic Steel[J]. Fatigue&Fracture of Engineering Materials & Structure,2000,23(15):747-758.

[15] JIANG Y Y , XU B,SEHITOGLU H . Three-dimensional Elastic-plastic Stress Analysis of Rolling Contact[J]. ASME Journal of Tribology,2002,124(4):699-708.

[16] SATOM, ANDERSON P M, ANDERSOND A. Rigney Rolling-sliding Behavior of Rail Steels[J]. Wear,1993,14(3): 159-172.