車輪不圓順輪軌滾動接觸行為的動態有限元仿真方法

（西南交通大學牽引動力國家重點實驗室，四川成都 610031）

1 研究背景

作為鐵路車輛中最為核心的走行部件——車輪，其系統的失效將會造成災難性后果及巨大的經濟損失。隨著列車速度的不斷提升和軸重的增加，車輪不圓順傷損問題日益嚴重，嚴重影響軌道-車輛零部件的服役壽命和車輛的乘坐舒適性，甚至威脅行車安全，已成為軌道交通領域備受關注和亟待解決的難點問題之一[1-3]。

目前，針對車輪不圓順條件下的輪軌關系研究通常采用多體動力學方法和有限元方法[4-6]。多體動力學方法大多基于Hertz接觸理論和穩態滾動假設，不能很好地描述輪軌系統的幾何、材料和接觸非線性，無法考慮輪軌相互作用的動態效應，且不能直接求解輪軌系統的應力/應變狀態。有研究表明，當研究對象無明顯結構變形（即頻率低于第一振型）時，多體動力學方法才是準確和有效的；而在頻率高于第一振型（大于84 Hz）時，其求解精度會顯著降低[7-8]。基于連續介質力學的有限元方法不僅可以考慮到局部變形，還可以考慮到強非線性、應變率效應和應力波傳播效應，被認為能夠更好地求解輪軌動態力學響應，近年來在求解輪軌滾動接觸行為研究中得到了越來越多的應用[9-12]。

然而，在采用有限元方法求解車輪不圓順引起的輪軌動態響應問題時，通常存在以下2個問題：①車輪踏面缺陷具有隨機性（包括尺寸、形狀、位置及分布的隨機性），造成幾何模型復雜，有限元建模難度很大；②車輪踏面缺陷處的網格精細劃分會引起有限元模型單元數量較大，造成計算時間過長，并且缺陷處不規則網格的畸變極易導致計算不穩定。針對上述問題，本文將車輪不圓順轉換為輪軌接觸位移不平順，引入偽隨機數，生成滿足多種隨機分布特性的輪軌接觸位移不平順激勵，并將其施加到不含缺陷的三維輪軌滾動接觸有限元模型中，極大地簡化了不圓順車輪-鋼軌間的動態接觸行為仿真研究。

2 車輪不圓順輪軌接觸數學模型

車輪不圓順（如車輪扁疤、踏面剝離、車輪多邊形等）是引起輪軌接觸不平順的主要原因之一。本文建立不圓順車輪和鋼軌滾動接觸數學模型，基于輪軌滾動接觸幾何關系，通過推導由車輪踏面缺陷引起的輪軌接觸不平順數學表達式，實現車輪不圓順向輪軌接觸位移不平順的轉換。

2.1 車輪扁疤

含扁疤車輪-鋼軌間滾動接觸情況如圖1所示。在低速運行情形下（圖1a），車輪圍繞扁疤前角轉動至扁疤與鋼軌平行后，緊接著圍繞扁疤后角開始轉動，直到越過整個車輪扁疤區域；在高速運行情形下（圖1b），車輪與鋼軌在扁疤區域內將會出現短暫的接觸損失（車輪脫離鋼軌），并在重力作用下以某一水平初速度做平拋運動，隨后扁疤后角撞擊鋼軌或直接越過扁疤繼續滾動。

本文采用具有最小曲率半徑的車輪輪廓運動軌跡描述實際的輪軌接觸位移不平順[13]。對于新扁疤，低速運行情形下的輪軌接觸位移不平順Zf（ξ）方程可表示為[13]：

而在高速運行情形下，結合平拋運動以及新扁疤車輪輪廓運動軌跡，其輪軌接觸位移不平順Zf（ξ）為：

式（1）～式（2）中，df，lf分別為新扁疤的深度和長度，其中df=lf2/ 8r；pw為車輪自重；pa為軸重；vl為水平初速度；g為重力加速度；ξ0為平拋運動結束后車輪與鋼軌接觸時新扁疤位置到扁疤起點的水平距離。

對于舊扁疤，低速運行情形下的輪軌接觸位移不平順Zr（ξ）可表示為[13]：

而高速運行情形下的Zr（ξ）則可表示為：

式（3）～式（4）中，ξ1，ξ3分別為舊扁疤前、后2個端點與中心區域的交點；a0為舊扁疤后端與中心區域的交點；a為確定過渡范圍的常量；df，lf分別為舊扁疤的深度和長度，其中df=alf2/4（2ar+1）。

2.2 踏面剝離

對于含踏面剝離的車輪，由于鋼軌頂面曲率的存在使得車輪剝離區域與鋼軌的接觸狀態變得較為復雜，因此分析中常將其簡化為規則的半橢球形。由于含踏面剝離車輪-鋼軌間動態接觸過程與扁疤類似，故可采用類似的分析方法推導規則踏面剝離造成的輪軌接觸位移不平順。踏面剝離的輪廓線可以簡化成一個圓柱體（車輪）與半橢球體（剝離）的交線，如圖2所示。此空間曲線可表示為：

式（5）中，（x，y，z）為車輪和踏面剝離在空間坐標系中的位置；ds，ls分別為踏面剝離的深度和長度。

圖1 含扁疤車輪-鋼軌間滾動接觸情況

圖2 含踏面剝離車輪與鋼軌接觸示意圖

通過幾何關系推導，可得到低速運行情形下踏面剝離引起的輪軌接觸不平順Z（y）表達式為：

而高速運行情形下的Z（y）表達式為：

式（6）～式（7）中，R1，R2分別為車輪名義滾動圓半徑和軌面曲率半徑；ws為踏面剝離寬度；m為輪對自重的一半；v為列車恒定速度；FS為作用在輪對上等效車體載荷的一半；θ為軌底坡角度；p為平拋運動階段的水平距離。

2.3 車輪多邊形

基于不圓順車輪半徑沿圓周方向的變化情況分為周期性和非周期性2種類型，其中非周期性多邊形車輪的徑向偏差可以看成不同諧波的疊加。周期性不圓順（三階）車輪-鋼軌間的滾動接觸過程如圖3所示。

對于周期性多邊形車輪，其徑向偏差?r可表示為[14]：

式（8）中，θ為車輪轉過的角度；A，N分別為某階車輪多邊形所對應的幅值和階數。

其中車輪多邊形波長λp和階數N之間的關系可表示為：

非周期性車輪多邊形磨耗可認為是一系列諧波的疊加，其中每個諧波具體指任意階的多邊形磨耗。因此，非周期性多邊形車輪引起的輪軌接觸位移不平順Z（θ）可表示為[15]：

式（10）中，φN為相位角，且在0～2π 之間服從多種隨機分布；AN為某階多邊形對應的幅值；M為車輪不圓順總的階次。

3 車輪踏面缺陷隨機性表征

如上所述，車輪踏面缺陷在尺寸、形狀、位置及分布上具有隨機性，本文通過選取合適的分布函數實現車輪踏面缺陷的隨機化描述，進而得到由車輪踏面缺陷引起的輪軌接觸隨機不平順位移激勵。

3.1 隨機分布函數

隨機數按生成方式整體上可分為真、偽隨機數2類[16]。其中真隨機數在生成之前無法預測，一般基于物理方法，通過一系列復雜過程產生；偽隨機數大多根據給定算法由數學方法產生，從原理上是可以預測的，通過選擇合理的計算方法及參數，產生具有良好特性的隨機數。此外，偽隨機數還具有生成效率高、實現簡單等特性。因此，本文采用偽隨機數表征車輪踏面缺陷的尺寸、形狀、位置及分布隨機特性。

偽隨機數的產生算法有很多，本文從計算速度、生成周期、獨立性和均勻性等方面綜合考慮，選取性能較好的線性同余法[17]。線性同余法是借助同余運算生成區間[0，1]內服從均勻分布的偽隨機數，其具體計算公式如下[18]：

圖3 周期性不圓順（三階）車輪-鋼軌間的滾動接觸過程

式（11）中，{rn}為[0，1]上的隨機序列；xn，xn-1分別為此隨機序列中的第n和第n-1個隨機數；b，c和M分別為乘子、增量和模，均是系統給定的常量；初值x0為種子。

由于車輪踏面缺陷在形成過程中的復雜性和不確定性造成其在車輪圓周上的分布情況極具隨機性，因此輪軌接觸位移不平順激勵通常可用多種隨機分布函數描述。例如：①在理想情況下，車輪踏面圓周上出現缺陷的位置是隨機的，且在任意一個位置出現的概率相等，故可用均勻分布的隨機數描述；②當車輪踏面缺陷的隨機性有多個獨立隨機因素疊加時，可認為每個隨機變量都滿足正態分布；③當車輪在某個位置已經出現缺陷，此后車輪每滾動到此處都會產生較大的振動沖擊，并在一定范圍內沖擊力會隨著車輪的滾動而不斷衰減，在此處再次產生缺陷的概率最大，往后依次遞減，此時可認為近似符合指數分布。因此，可以通過對標準均勻分布、正態分布、指數分布等進行相應變換使其滿足特定區間內的隨機分布類型，實現車輪踏面缺陷的隨機化描述。

3.2 輪軌接觸隨機不平順位移激勵

基于踏面缺陷引起的車輪不圓順與輪軌接觸位移不平順之間的轉換關系，本文采用線性同余法生成滿足各種隨機分布函數的輪軌位移不平順激勵。限于篇幅，這里僅以踏面剝離為例，給出不同隨機度λz下單個踏面剝離形貌及相應的輪軌接觸位移不平順曲線，如圖4所示。

4 有限元仿真方法及驗證

4.1 有限元模型

借助Hypermesh軟件建立包含軌下基礎結構的三維輪軌滾動接觸有限元仿真模型，如圖5所示。車輪為S1002CN型踏面，半徑為430 mm；鋼軌為CN60型，軌底坡為1 :40。整個模型由彈簧-阻尼離散單元（扣件）和8節點實體單元2部分組成。為減小計算量，將彈簧阻尼單元連接區域和輪軌接觸區域的網格進行細化（4 mm×4 mm），其余網格進行適當過渡。整個模型包含1 177 793個節點和1 001 343個單元（包括實體單元和離散單元）。

圖4 不同隨機度λz 下單個踏面剝離的形貌與輪軌接觸位移不平順激勵

圖5 三維輪軌滾動接觸有限元仿真模型

車輪和鋼軌均采用*MAT_PLASTIC_KINEMATIC材料模型描述，該模型可以模擬材料的隨動強化特性，且可通過Cowper-Symonds本構關系考慮輪軌材料的應變率效應。車軸采用*MAT_ RIGID材料模型描述，軌下基礎結構各部分采用*MAT_ELASTIC材料模型描述，彈簧和阻尼系統采用*MAT_SPRING_ELASTIC和*MAT_DAMPER_VISCOUS材料模型描述。仿真中輪軌系統各部件的力學性能參數如表1所示。

為實現輪對運動過程中的穩定與平衡，仿真中對路基底部和車軸端部所有節點分別實施全約束和軸向平動約束，給輪對設置相應的初始平動和轉動速度實現其在鋼軌上的滾動。根據歐洲標準，在車軸兩端分別施加P1=77.56 kN和P2=110.41 kN 2個集中力等效軸重載荷（17 t）。車輪與鋼軌、車軸之間均定義基于罰函數法的Automatic_Surface_To_Surface，軌下基礎結構各部件之間定義Tied_Surface_To_Surface接觸。如上所述，為避免車輪踏面缺陷帶來的有限元直接建模諸多問題，仿真中將車輪踏面缺陷引起的輪軌接觸位移不平順激勵作為初始邊界條件，通過將關鍵字*PRESCRIBED_MOTION_SET輸入到鋼軌頂面的相應節點上，并采用強迫節點運動方式實現對車輪踏面缺陷的模擬。

表1 仿真中輪軌系統各部件的力學性能參數

有限元仿真中的積分算法通常包括隱式和顯式2種。顯式算法適用于求解波的傳播、非線性動力學等高頻、瞬態問題，其積分時間步一般很小，且求解速度很高，不存在收斂性問題。其瞬態動力學分析方程為：

式（12）中，[M]，[C]，[K]分別為結構質量、阻尼和剛度矩陣；{F（t）}為結構外載荷；為結構節點加速度矢量；為結構節點速度矢量；為結構節點位移矢量。

由于高速輪軌滾動接觸涉及高頻、瞬態動力學分析，且必須考慮動態效應的影響，故本文選用顯式積分算法。

4.2 方法應用及驗證

基于上述動態有限元仿真方法（為與基于不圓順車輪真實幾何輪廓的直接建模方法區別，下面稱為位移激勵法），以車輪扁疤為例，采用LS-DYNA軟件求解不圓順車輪-鋼軌間動態響應與接觸行為，并與直接建模得到的仿真結果相比較，以驗證該方法的計算準確性和精度。

在特定工況（v=200 km/h，l=40 mm）下，基于位移激勵法和直接建模法計算的輪軌垂向接觸力-時間曲線及差異性[19]如圖6所示?？梢钥闯觯逯盗ι杂胁町愅猓@2種方法的計算結果吻合較好；不同速度和扁疤長度下2種方法得到的峰值力相對誤差均不超過17.8%。如前所述，高速輪軌關系研究中還必須考慮輪軌動態效應，特別是輪軌系統材料的應變率效應。圖7所示為含車輪扁疤條件下，材料應變率效應對輪軌系統應力狀態及分布規律的影響。可以發現，在不考慮應變率效應時的輪軌接觸應力峰值分別為583.3 MPa和549.4 MPa，而在考慮應變率效應時的輪軌接觸應力峰值為648.4 MPa和626.8 MPa，相應增加了65.1 MPa和77.4 MPa。這表明應變率效應能夠顯著提高輪軌接觸von Mises應力，但是對其分布情況幾乎沒有影響。

圖6 車輪扁疤引起的輪軌垂向接觸力響應

圖7 車輪扁疤引起的輪軌接觸von Mises應力狀態及分布規律

上述含車輪扁疤的輪軌動態接觸行為仿真示例表明，本文提出的將車輪不圓順轉換為輪軌接觸位移不平順的動態有限元仿真方法是可行、有效的。在考慮輪軌系統強非線性特征和動態效應影響的同時，能夠較好地表征缺陷的隨機性特點，并且有效地求解輪軌接觸應力/應變狀態，以及分析高速不圓順車輪-鋼軌間的滾動接觸行為。

5 結語

針對采用有限元方法求解含踏面缺陷車輪-鋼軌間滾動接觸行為時的諸多困難，本文通過建立車輪不圓順條件下的輪軌接觸數學模型，提出了用輪軌接觸位移不平順代替車輪踏面缺陷的等效方法；基于線性同余法生成的偽隨機數，通過選取合適隨機函數實現了車輪踏面缺陷在分布位置、幾何形狀/尺寸和數量等方面的隨機化表征，以初始位移激勵形式施加到輪軌滾動接觸仿真模型中，在考慮輪軌強非線性特征和動態效應影響的同時，有效地求解輪軌動態力學行為（應力/應變狀態等）；最后，通過與直接建模得到的仿真結果相比較，得出位移激勵法與直接建模法的輪軌間動態響應吻合較好，驗證了基于位移激勵的車輪踏面缺陷有限元模擬方法是可行、有效的。該仿真方法為高速輪軌滾動接觸關系的研究提供了一種新的思路。