周新建,盧 勇,王成國,肖 乾,,張 懷
(1 華東交通大學 載運工具與裝備教育部重點實驗室,江西南昌330013;2 中國鐵道科學研究院 鐵道科學技術研究發(fā)展中心,北京100081;3 中國科學院 研究生院 計算地球動力學實驗室,北京100049)
目前國產(chǎn)大功率機車已批量下線并投入運用,促進了我國民族工業(yè)發(fā)展,形成了以重點企業(yè)為核心、配套企業(yè)為骨干、輻射上百家相關企業(yè)的國內機輛裝備設計制造體系。大功率機車走行部,尤其是車輪這類關鍵零部件。實現(xiàn)國產(chǎn)化的最重要步驟就是需要對目前大功率機車的現(xiàn)場運行情況有非常清晰地了解,掌握全面的基礎數(shù)據(jù)。輪軌接觸應力分析是輪軌接觸問題的基礎。對于機車車輪來說,必須明確其在不同工況下與鋼軌的接觸特性,因此,完成大功率機車的輪軌接觸應力計算分析就迫在眉睫。
接觸問題通常采用增量方法求解,與一般的非線性問題區(qū)別在于每一個增量步中都必須將接觸的約束條件作為定解條件和校核條件判斷接觸狀態(tài)。也就是說,計算時,首先得到和時間t+Δt位形內平衡條件相等效的虛位移方程,然后在此基礎上求得其T.L.或U.L.格式方程,一般式描述如下:
式中M是節(jié)點質量矩陣,K為節(jié)點剛度矩陣,P為節(jié)點載荷矩陣,μ為節(jié)點位移,M取零值為靜態(tài)接觸問題。這里需要特別指出的是根據(jù)拉格朗日乘子法和罰函數(shù)法引入約束條件的不一樣,接觸問題有限元求解方程的具體形式也有區(qū)別。對于這些非線性有限元方程組,可以采用顯式解法或隱式解法進行求解。顯式解法求解接觸問題時,通常采用二步形式的中心差分法進行遞推[1],即有:
而隱式解法常常采取Newton-Raphson方法進行迭代求解,N-R的迭代公式為[2]:
顯式解法不僅可以避免每一增量步的迭代步驟,而且可以避免摩擦滑動接觸狀態(tài)時非對稱剛度矩陣的求逆運算,但在每個時間步長內沒有誤差檢查和迭代步驟,計算誤差不能控制。而隱式解法的優(yōu)缺點與顯式求解正好相反,能保證求解精度但計算時間相對較長。為了保證良好的計算精度,輪軌接觸有限元計算采用隱式求解法。隱式算法的計算速度主要由計算機內存決定,采用并行計算能完成單機無法計算的大規(guī)模有限元模型,考慮到計算模型規(guī)模比較大,用單機無法勝任,所以采用有限元并行計算。
1.2.1 輪對有限元模型
大功率機車車輪幾何模型是以HXD2機車車輪參數(shù)繪制,直徑為1 250 mm,采用磨耗型(JM)踏面。由于網(wǎng)格細化是保證計算精度的主要手段之一[3],文獻[4]以對有限元兩次分析的結果進行比較,如果二者相差小于2%,說明當前的網(wǎng)格密度和單元類型足以保證工程分析精度,指出在輪軌計算中,單元尺寸須小于等于0.75 mm。為了有較好的精度,對車輪可能接觸的區(qū)域進行局部網(wǎng)格細化,細化區(qū)網(wǎng)格尺寸為0.5 mm,采用一階六面體減縮積分單元(C3D8R)進行離散,整個輪對產(chǎn)生733 995個單元和758 126個節(jié)點。
1.2.2 鋼軌有限元模型
鋼軌采用的是75 kg/m國內標準軌,縱向長度取5 m。根據(jù)輪對細化相同的原理,對鋼軌可能的接觸區(qū)域進行網(wǎng)格細化,細化區(qū)最小單元尺寸為0.5 mm。鋼軌也采用一階六面體減縮積分單元(C3D8R)進行離散,兩段鋼軌共產(chǎn)生252 030個單元、267 430個節(jié)點。
1.2.3 模型的材料、邊界條件以及接觸面的定義
輪軌材料模型為經(jīng)典的雙線性彈塑性本構模型,材料本構曲線如圖1。軌距為1 435 mm,軌底坡為1:40,軸重Q為25 250 kg,重力加速度為9.81 m/s2,泊松比υ為0.3。根據(jù)UIC 510-5[5]標準規(guī)定的計算載荷,對應直線運行時的受力狀態(tài)。車輪受到垂直方向力的作用:
式中Q是單個車輪承受的質量,g為重力加速度。通過公式(4)可求得單個軸箱上所要承受的載荷為154 814 N。加載時以集中應力形式加載在輪對兩端軸箱位置處,同時約束鋼軌底面的所有自由度。接觸面定義為柔性體間的面面接觸,接觸算法選用擴展的拉格朗日算法,摩擦模型選用罰摩擦模型。輪軌接觸有限元模型見圖2。由于輪軌接觸有限元模型規(guī)模較大,采用中國科學院研究生院計算地球動力學實驗室的網(wǎng)絡集群并行計算環(huán)境,該并行系統(tǒng)計算節(jié)點的配置為2顆4核E5520處理器、8G DDR3內存和146GSAS接口的HP Proliant BL280c G6刀片服務器,根據(jù)ABAQUS6.9.1并行求解器的特點,采用12個節(jié)點共96核進行接觸有限元并行計算。
圖1 輪軌非線性本構曲線
圖2 輪軌接觸有限元網(wǎng)格圖
在不考慮超高和沖角的條件下,對車輪橫移量為-9,-6,-3,0,3,6,9 mm 7種工況下的輪軌接觸情況進行計算,獲得了不同橫移情況下的接觸斑形狀、接觸斑面積、車輪的表面最大接觸壓應力、車輪最大Mises應力以及車輪的最大等效塑性變形量等。
表1所示的是大功率機車車輪橫移為0時的輪軌接觸計算結果。
表1 橫移為0時輪軌接觸計算結果
圖3為輪對橫移為0 mm時的輪軌接觸Mises應力分布云圖,車輪最大Mises應力出現(xiàn)在接觸表面以上2.5 mm處的車輪踏面內。圖4給出了輪軌最大等效塑性變形分布云圖,從圖3和圖4比較可以發(fā)現(xiàn)鋼軌的Mises應力大于車輪的,但是由于車輪的屈服應力小于鋼軌的(如圖1所示),所以在車輪踏面出現(xiàn)了較大的塑性變形區(qū),最大塑性變形點出現(xiàn)位置與車輪最大Mises應力出現(xiàn)位置一致,車輪接觸表面也有塑性變形,車輪塑性變形區(qū)最大縱向長度為16.5 mm,最大橫向長度為6.5 mm,最大垂向長度為5.7 mm。從圖5可以看到接觸斑橫向長度向背離輪緣側延伸,與Hertz理論的橢圓形接觸斑有較大差異。
圖3 橫移為0時輪軌接觸的彈塑性Mises應力分布云圖
圖4 橫移為0時等效塑性 變形分布
圖5 橫移量為0時 接觸斑形狀
參照文獻[6],提取不同橫移時大功率機車的輪軌接觸計算結果(如表2所示)。車輪踏面和軌頭都是由多段不同半徑的圓弧構成,表 2中所列的 Rwheelprof和Rrailprof分別是指接觸斑所處的車輪踏面和鋼軌軌頭的圓弧半徑,從表2中所列的不同橫移時的接觸斑形狀可以看到,輪軌接觸具有明顯的非Hertz接觸特點。
表2 不同橫移計算結果數(shù)據(jù)表
從圖6可以看出接觸斑橫向長度和縱向長度的變化呈現(xiàn)互補特性,即橫向長度增大時縱向長度變短,反之亦然。由圖6和圖7可以看到,在橫移量為-3~-6 mm時輪軌接觸面積比較大,當橫移量絕對值超過6 mm時,輪軌接觸面積隨橫移增大而變小,而且接觸面積隨橫移量的變化曲線與接觸斑橫向長度隨橫移量的變化曲線有相同的變化規(guī)律。
從圖8可以發(fā)現(xiàn),大功率機車輪對橫移量的絕對值在6 mm以內變化時,Mises應力變化不大,當橫移量從6 mm增加到9 mm時,車輪最大Mises應力迅速升高。這是由于當橫移量為9 mm時,車輪踏面接觸圓弧最小曲率半徑與鋼軌軌頭接觸圓弧最小曲率半徑分別為14 mm和15 mm,接觸圓弧的變小使得接觸壓力增大,而且,此時輪軌接觸區(qū)域為一小半徑傾斜曲面,曲面上的點在車輪徑向上線速度相差較大,這會導致車輪和鋼軌有相對滑動產(chǎn)生,加劇車輪的磨損;接觸應力迅速升高也加速了車輪側磨的發(fā)生,所以應減少橫移量超過6 mm的時候;最大接觸壓力隨橫移量的變化幅度相比最大Mises應力的變化幅度大,橫移為6 mm時接觸壓力最小。
從圖9可見,輪對在不同橫移時,車輪踏面都發(fā)生局部塑性變形。在橫移為6 mm時最小,橫移量從6 mm變化到9 mm時,最大等效塑性變形量迅速變大。
圖6 車輪接觸斑長度隨橫移量的變化
圖7 車輪接觸斑面積隨橫移量的變化
圖8 車輪踏面體內外接觸應力隨橫移的變化
圖9 車輪最大塑性變形隨橫移量的變化
(1)車輪踏面內出現(xiàn)塑性變形,塑性變形從車輪踏面內延伸至接觸表面。對橫移為0時車輪踏面的分析可知車輪塑性變形區(qū)最大縱向長度為16.5 mm,最大橫向長度為6.5 mm,最大垂向長度為5.7 mm,而且不同橫移時車輪踏面都出現(xiàn)塑性變形,可見車輪踏面塑性變形不可避免。
(2)隨著車輪橫移量的不同,橫向接觸斑長度與接觸面積有著相同的變化規(guī)律,多數(shù)情況下的輪軌接觸斑形態(tài)與Hertz理論的橢圓假設有較大差別。
(3)大功率機車輪對橫移量的絕對值在6 mm以內變化時,Mises應力變化不大;當橫移量從6 mm增加到9 mm時,車輪最大Mises應力迅速升高,此時車輪和鋼軌接觸區(qū)的最小圓弧曲率半徑分別為15 mm和14 mm,小曲率半徑圓弧相接觸是造成接觸應力升高的主要原因。應盡量減少機車橫移量超過6 mm的情況。
(4)運用接觸有限元并行計算技術可以解決輪軌接觸有限元計算的計算規(guī)模和計算精度的矛盾,這也是未來用有限元計算輪軌接觸狀態(tài)的趨勢。輪軌接觸應力應變數(shù)據(jù)是研究大功率機車車輪的基本依據(jù),計算結果為大功率機車車輪國產(chǎn)化提供理論參考。
[1] TED B,L IU Wingkam,BRIAN M.Nonlinear Finite Element s for Continua and Structures[M].New York:John Wiley&Sons,2000.
[2] 羅喜恒,肖汝誠,項海帆.幾何非線性問題求解的改進算法[J].公路交通科技,2005,22(12):75-77.
[3] 石亦平,周玉蓉.ABAQUS有限元分析實例詳解[M].北京:機械工業(yè)出版社,2008.
[4] XiaoQian Wang Chenguo GuoGe.”The Research of Parallel Computing for Large-scale Finite Element Model OF Wheel/Rail Rolling Contact”.Proceedings of 2010 3rd IEEE International conference on Computer Science and Information Technology,Chengdu,pp:254-257,July.2010.
[5] UIC 510-5.鍛造車輪和軋制車輪[S].國際鐵路聯(lián)盟.2003.
[6] Coenraad Esveld.Modern Railway Track(Second Edition)[M].The Netherlands:MRT-Productions,2001.