三種非赫茲滾動接觸模型的對比研究

王 平，周佳儀，王攀杰，陳 爽，安博洋

(1.西南交通大學 高速鐵路線路工程教育部重點實驗室，四川 成都 610031；2.西南交通大學 土木工程學院，四川 成都 610031；3.成都工業職業技術學院 軌道交通學院，四川 成都 610031)

鐵路運營速度的不斷提高使得輪軌服役環境日益嚴苛，鋼軌波浪形磨損、鋼軌剝離、車輪多邊形磨耗等輪軌傷損現象也更為嚴重[1]。因此，更為準確地模擬輪軌滾動接觸力學行為，對于預測輪軌滾動接觸疲勞和磨耗具有十分重要的意義。以往的模擬研究中[2-3]，為了提高計算效率，常采用基于橢圓假設的赫茲接觸理論。然而，實際情況下的輪軌接觸斑常為非橢圓形狀，使用赫茲理論會造成較大的計算誤差。輪軌非赫茲滾動接觸模型可以更準確地描述輪軌的接觸斑形狀及接觸應力分布，從而能夠更好地預測輪軌磨耗及滾動接觸疲勞。

在模擬輪軌非赫茲滾動接觸的研究中，業界學者開發了兩類準確的計算模型。其一是通過有限元方法建立的數值模型，但是此種方法需要對輪軌結構進行精細的網格劃分，因此在計算過程中需要進行龐大的矩陣計算，極大地限制了運算速度[4-5]。另外一種是基于半空間方法建立的邊界元方法，其中最為典型的代表是Kalker開創的三維彈性體非赫茲滾動接觸理論及其數值模型CONTACT[6-7]。該理論是目前業界最廣為接受的理論，但是其計算效率因不斷迭代而難以提高，因此尚未廣泛應用于輪軌損傷的計算。

為了更高效地計算非赫茲滾動接觸問題，學者們發展了一系列快速解法。Piotrowski等[8]提出了虛擬滲透的概念，通過引入衰減系數以忽略材料彈性變形的影響，提出了一種快速求解輪軌法向接觸問題的近似方法(下文簡稱為KP模型)。LINDER等[9]也提出了類似的思路，但通過不同方法來計算輪軌接觸斑和接觸應力。Liu等[10]考慮搖頭角的影響對KP模型進行了更深一步的擴展。文獻[11]同樣在KP模型的基礎上發展了一種考慮彈性變形與法向壓力分布關系的改進模型，使計算不再依賴于接觸原點。另一種基于虛擬滲透方法的經典模型是文獻[12]提出的半赫茲接觸模型，該模型將滲透區域沿滾動方向劃分為數個條帶，并假定每個條帶滿足赫茲接觸條件，從而求解出接觸區域以及應力分布，并發展了STRIPES的數值程序。與STRIPES類似，文獻[13]考慮材料彈性變形的影響，提出了一種新的方法——ANALYN來解決輪軌法向接觸問題[13]。而以上簡化模型切向問題的處理均采用改進的FASTSIM方法[8,12]或FaStrip模型[13]。其中改進FASTSIM方法的核心是將非橢圓接觸斑等效為一個或多個局部橢圓，以求解出新的柔度系數替代原方法中根據橢圓假設計算的系數。

針對上述非赫茲接觸模型，文獻[14]比較了KP、LINDER和STRIPES三種方法求解S1002-UIC60接觸工況時法向與切向結果的差異。文獻[15]詳細分析了KP、STRIPES、LINDER三種模型與FASTSIM程序在曲線地段進行切向計算的差異。文獻[16]也對典型接觸模型在法向以及切向方面的計算差異進行了對比，并分析了不同接觸模型對輪軌傷損模擬的影響。文獻[17]針對道岔轉轍器區以CONTACT程序為參照，對比分析了不同赫茲接觸模型，包括KP模型、STRIPES模型以及ANALYN模型在法向接觸行為以及切向接觸行為方面的計算差異。上述研究均表明，簡化的非赫茲接觸模型的精度對輪軌型面的依賴性很大，且只針對有限的接觸工況進行對比，故其適用范圍存在一定的局限性。因此，仍然需要一個系統、全面的對比研究，考慮不同型面類型，對大量工況進行計算，比較不同的簡化非赫茲模型在不同接觸工況下的計算差異。

本文選取我國高速鐵路常見的兩種標準車輪型面(LMA、S1002CN)與標準CHN60鋼軌匹配時的兩種典型工況進行詳細接觸計算，并以CONTACT程序的結果作為參照。借助UM軟件輸出大量工況，采用統計的思想，系統全面地調查不同模型在計算輪軌滾動接觸行為時的精度與適用性，為輪軌損傷計算的模型選擇提供可靠的依據。

1 三種非赫茲接觸簡化模型描述

1.1 KIK-PIOTROWSKI計算模型

KIK與PIOTROWSKI提出了一種快速、近似解決輪軌法向接觸的新方法(以下簡稱“KP模型”)。在輪軌接觸中，以接觸點為坐標原點定義局部坐標系，x為縱向，正方向代表車輪滾動前進方向，y為橫向，可以通過右手法則確定。該方法在法向接觸建模中，假設兩物體之間可以相互穿透，采用經驗系數(ε=0.55)與滲透量的乘積來確定接觸區域，并認為接觸區域內的壓應力分布應與縱向接觸邊界一致。其縱向接觸邊界a(y)及法向接觸應力分布p(x,y)為

( 1 )

( 2 )

式中：f(y)為接觸點附近沿橫向的法向間隙；δ為滲透量；R為接觸點處車輪半徑；p0為接觸原點的法向接觸應力。

切向接觸問題的求解則采用修正的FASTSIM模型進行處理。該方法的核心思想是將非橢圓區域等效為一個橢圓，根據該橢圓的具體參數確定FASTSIM中的柔度系數。詳細介紹可參見文獻[8]。

1.2 STRIPES計算模型

同樣基于輪軌虛擬穿透的基本原理，AYASSE和CHOLLET在一系列推導過程的基礎上發展了半赫茲的簡化模型STRIPES。將接觸斑沿滾動方向分解為若干條帶，并假設每個條帶的特征與赫茲接觸理論相同。通過修正接觸點區域的相對曲率，求解輪軌滲透量確定接觸區域內法向接觸應力的分布。每個條帶中的縱向半軸長度a(y)以及接觸斑內法向接觸應力p(x,y)為

( 3 )

( 4 )

切向接觸計算時仍采用修正的FASTSIM模型，不同于KP模型將整個接觸斑等效為一個橢圓，STRIPES模型視每個條帶均通過一個虛擬橢圓中心，在赫茲理論中使用條帶中心處的相對曲率來確定局部虛擬橢圓長短半軸，由局部虛擬橢圓參數計算每個相應條帶的柔度系數，隨后根據原始FASTSIM方法進行切向接觸應力計算。具體詳見文獻[12]。

由于運動副連續接觸，使運動副的銷軸與套圈中心始終保持運動副間隙的距離.可以將各運動副間隙假設為無質量的剛性桿.剛性桿長度為運動副間隙量，即r=(r1,r2,r3,r4)，各剛性桿與X軸正方向的夾角稱為運動副間隙角，記為α=(α1,α2,α3,α4).用剛性桿表示運動副間隙，運動副間隙角反映銷軸與套圈的接觸位置的變化.

1.3 ANALYN計算模型

SICHANI提出了一種不同于虛擬滲透法的近似分析方法，該方法不再忽略表面變形的影響，而是通過接觸面之間的間隙估計接觸表面變形，并通過修正滲透量求得接觸區域。該方法認為接觸區域的邊界應滿足：接觸物體之間的法向間隙等于物體之間法向剛體穿透量與接觸表面法向彈性變形量之差。接觸斑的縱向和橫向邊界及法向接觸應力為

( 5 )

( 6 )

( 7 )

式中:α(y)、β(y)、n(y)、r(y)為每個條帶處由赫茲理論計算所得的參數；A(y)、B(y)為每根條帶上的赫茲計算參數，即條帶中心處的曲率；p0(y)為每個條帶上的最大接觸壓力。

切向接觸問題的解決以KALKER的條帶理論和線性理論為基礎，發展成一種名為FaStrip的新模型。該方法在計算過程中使用統一的柔度系數，根據每個條帶處的相對曲率參數確定虛擬橢圓的長短半軸，并判斷每個條帶處于黏著區或滑動區，從而獲得對接觸滑動區域分布的準確估計。詳細介紹參見文獻[13]。

2 典型工況下輪軌非赫茲滾動接觸分析

2.1 典型工況的選取

為了對比上述三種非赫茲快速模型的計算差異性，以KALKER的CONTACT程序的計算結果作為標準參考。所考慮的典型工況為LMA、S1002CN兩種車輪型面分別與CHN60鋼軌匹配時輪軌接觸斑在鋼軌頂面的分布情況。其中，輪對橫移量變化范圍為-3～6 mm，LMA和S1002CN車輪名義滾動圓半徑分別為430、460 mm，輪對內側距和軌距分別為1 353、1 435 mm，軌底坡設置為1/40，輪對軸重14 t。特別指出，由于本文所考慮的簡化模型均假設輪軌接觸斑沿縱向滿足赫茲型接觸假設，故本節的模擬工作忽略了輪對搖頭角的影響，而輪對搖頭角引起的橫向蠕滑率則在下文中加以考慮。

不同輪軌型面匹配時的接觸斑形狀以及法向應力分布隨橫移量的變化見圖1。為了便于比較分析，選取LMA型面非赫茲程度較小、接觸斑形狀接近橢圓，即橫移量為0 mm時作為非極端工況代表。由圖1(j)可知，S1002CN型面在6 mm橫移量下接觸斑非赫茲程度比較大，因此選取該工況作為極端工況代表。從接觸斑形狀、法向接觸應力分布以及切向黏滑區域分布等方面對兩種工況進行計算，分析不同非赫茲程度下三種模型的計算精度以及適用性。

圖1 不同輪軌型面匹配下，輪軌接觸斑形狀以及法向應力分布隨輪對橫移量(-3～6 mm)的變化

2.2 輪軌法向非赫茲接觸分析

圖2為不同接觸工況下三種簡化接觸模型與CONTACT程序計算的接觸斑面積、最大法向接觸應力以及接觸區域內橫向曲率分布，其中，圖2(a)～2(c)為LAM型橫移量為0 mm時的接觸工況，圖2(d)～2(f)型為S1002CN橫移量為6 mm時的接觸工況。如圖2(a)～2(c)接觸斑形狀非赫茲程度較小時，三種簡化接觸模型得到的接觸斑形狀與CONTACT計算結果基本吻合。而隨著非赫茲程度的增大，快速模型的計算誤差也逐漸增大，但ANALYN的計算結果始終保持與CONTACT程序最為接近。由圖2可以看出，KP模型總是高估了接觸斑面積，致使相應法向接觸應力偏小，說明以0.55倍的虛擬滲透量值作為實際滲透量并不完全準確。根據KP模型壓力分布公式計算出的壓力分布形狀類似于接觸斑的邊界，在接觸斑非赫茲程度較大時會產生較大的誤差。在給定的載荷下，接觸斑面積越大，接觸斑內分布法向應力越小，因此KP計算所得分布法向應力總小于CONTACT。此外，值得注意的是，LMA型面的橫向相對曲率分布比較平滑，而S1002CN廓形的橫向相對曲率分布存在較大的跳躍波動，甚至會出現負值，違背了赫茲理論中曲率恒定的假設，進而對依賴于局部曲率進行計算的STRIPES模型、ANALYN模型產生一定影響，也是造成接觸斑非赫茲程度較大的主要原因。

圖2 不同接觸模型計算不同工況時法向接觸對比

2.3 輪軌切向非赫茲接觸分析

輪軌接觸斑上相互接觸的質點對之間存在的相對滑動會引起接觸區域內的切向力，又叫做蠕滑力，是引起輪軌磨耗接觸疲勞和傷損的主要因素，也是車輛軌道動力學研究的重要輸出量。因此，需要對不同型面接觸產生的切向力以及黏著區域分布進行分析。分析過程中忽略搖頭角以及側滾角等的影響，因此橫向蠕滑力為零，縱向蠕滑率ξxL,R、自旋蠕滑率ξnL,R為[18]

( 8 )

( 9 )

式中：rL,R為左右輪瞬時滾動圓半徑；r0為名義滾動圓半徑；δL,R為左右輪接觸處的接觸角。圖3、圖4是兩種典型工況下蠕滑黏著區域的分布，接觸斑內顏色的深淺代表蠕滑應力的大小，零點代表接觸點的位置。每種接觸模型計算所得切向力數值見表1。

圖3 LMA型面橫移量0 mm時輪軌黏滑分布

圖4 S1002CN型面橫移量6 mm時輪軌黏滑分布

表1 不同型面不同橫移量下的縱向、橫向以及總蠕滑力

3 基于統計思想的輪軌非赫茲接觸簡化模型適用性分析

為了檢驗以上簡化模型對于車輛動力學工況的適用性，借助適用于多體系統運動學、動力學建模和仿真計算的UM軟件構建CRH2單節車輛模型進行仿真分析，從而獲得接觸計算所需的基礎參數。圖5為單節車輛模型。

圖5 單節車輛模型

根據多體動力學理論，將CRH2型車視為非線性多剛體系統，單節車輛由15個剛體組成。構架與輪對通過一系懸掛裝置連接，車體與構架通過二系懸掛裝置連接，采用不同的鉸接約束車輛中各剛體的運動。采用LMA、S1002CN兩種我國常用廓形車輛與CHN60鋼軌匹配，設置速度為250 km/h，分別通過不平順直線以及平順曲線地段，并設置兩種不同摩擦系數(f=0.3、f=0.5)，從UM軟件中輸出蠕滑率等計算參數，將其作為輸入參數導入CONTACT程序以及MATLAB編制的KP模型、STRIPES模型、ANALYN模型中進行計算，將不同模型的計算結果與CONTACT進行對比，采用統計中累計分布的思想對誤差進行處理，對比簡化接觸模型的計算精度，計算流程見圖6。需要強調地是，本節的研究目的是通過動力學計算獲得大量的輪軌動態響應，以便于從統計的角度分析不同非赫茲接觸模型的適用性，因而僅考慮了一種車型的參數。

圖6 計算流程

3.1 直線地段

模擬車輛通過不平順直線軌道，研究車輛不平順對不同接觸模型計算精度的影響。在直線軌道上施加德國高速不平順譜，獲得S1002CN、LMA車輪型面與CHN60鋼軌匹配時的橫向、縱向、自旋蠕滑率以及橫移量等參數，其概率分布見圖7。然后將已獲得的接觸參數分別代入三種簡化接觸模型以及CONTACT程序進行法向及切向接觸計算。

圖7 UM軟件輸出不同型面接觸參數概率分布

圖8 不同接觸模型計算LMA型面接觸時累積誤差分布

使用UM軟件輸出不同接觸型面匹配時的接觸參數作為接觸模型的輸入量，采用三種非赫茲接觸模型進行計算，并以CONTACT模型計算結果作為參考值，得到大量工況下幾種接觸參數的計算結果，采用統計的思想求得三種模型對CONTACT的誤差累積概率分布。橫坐標表示與CONTACT模型計算結果相比產生誤差的數值，縱坐標表示達到該誤差值工況的出現概率。UM軟件輸出不同型面接觸參數概率分布見圖7。由圖7可知，車輛在不平順直線地段行駛時所得蠕滑分布較為集中，且S1002CN型面產生的自旋值會略大于LMA型面。根據2節的結果可知，由于型面接觸范圍內曲率存在較大的波動，S1002CN型面在大蠕滑條件下會產生較大誤差。采用不同接觸模型計算LMA、S1002CN型面接觸時蠕滑力累積分布分別見圖8、圖9。由圖8、圖9可知，S1002CN廓形下計算得到的各項接觸參數累積誤差要大于LMA廓形。在不同接觸模型對同一接觸參數的計算中，圖9表明ANALYN模型總能獲得最優結果；KP模型雖然預測詳細接觸解時存在一定的不準確性，但在計算接觸參數時多數情況下可以得到不錯的結果；而STRIPES模型相對最差。此外摩擦系數的變化不會影響接觸參數以及法向接觸計算，因此接觸斑面積、最大法向接觸應力以及總法向力不變。而由圖8、圖9可以看出，提高摩擦系數對切向計算的影響不大。

圖9 不同接觸模型計算S1002CN型面接觸時蠕滑力累積誤差分布

3.2 曲線地段

模擬車輛在50 m的直線軌道上行駛，然后通過400 m過渡曲線，進入半徑為3 000 m的曲線。通過UM軟件獲得S1002CN、LMA車輪型面與CHN60鋼軌匹配時的橫向、縱向、自旋蠕滑率以及橫移量等參數，其概率分布見圖10。將已獲得的接觸參數分別代入三種簡化接觸模型以及CONTACT程序進行法向及切向接觸計算。

圖10 UM軟件輸出不同型面接觸參數概率分布

由圖10可見，車輛由直線逐漸過渡到曲線地段時蠕滑分布較為分散。不同接觸模型計算LMA型面、S1002CN型面的接觸參數累積誤差分布分別見圖11、圖12。由圖11、圖12可知，S1002CN型面接觸參數的計算誤差要大于LMA型面，且通過對比圖8、圖9可以看出，曲線地段的接觸計算誤差略大于不平順直線地段，原因在于曲線地段會產生較高比例的高蠕滑。與直線地段類似，隨著摩擦系數的提高，各模型切向計算的精度差異不大。在曲線地段，KP模型的計算精度與ANALYN模型基本一致，甚至高于ANALYN。

圖11 不同接觸模型計算LMA型面的接觸參數累積誤差分布

圖12 不同接觸模型計算S1002CN型面的接觸參數累積誤差分布

以上內容均從計算精度的角度對不同非赫茲滾動接觸模型進行了對比分析，而車輛動力學求解和輪軌傷損計算要求接觸模型的求解效率越快越好，不同接觸模型在同一臺筆記本電腦上計算同一工況所需時間見表2。其中，CONTACT程序所需時間最長，約為119.2 s；三種快速計算模型中KP模型的計算時間最短，約為CONTACT程序的1/397；ANALYN模型與STRIPES模型所需時間較KP模型較長。需要指出地是，本文所用的 KP模型、STRIPES模型以及ANALYN模型都是編寫的MATLAB程序。而CONTACT是在Kalker原始算法基礎上經過Vollebregt優化后的求解程序[19]。

表2 不同模型同一工況計算時間對比

4 結論

本文采用目前最常見的三種非赫茲接觸模型從輪軌法向接觸和切向接觸兩個方面對我國高速車輪廓形LMA以及S1002CN與CHN60鋼軌匹配時進行計算分析，研究了不同非赫茲接觸模型計算差異，根據計算結果得到以下結論:

(1)三種簡化模型的計算精度在處理LMA-CHN60接觸工況時最優，該接觸工況下，接觸斑形狀接近于橢圓形，非赫茲程度較小，三種快速計算模型得到的接觸斑形狀、接觸區域面積以及最大法向應力值與CONTACT計算結果基本相同。

(2)S1002CN廓形存在明顯的局部曲率突變，更易產生非赫茲程度較大的接觸斑，簡化模型的計算精度在處理此類接觸問題時較差。ANALYN模型在計算黏滑區域分布方面與CONTACT的結果最接近，因此計算詳細局部黏滑分布時建議采用ANALYN模型。

(3)宏觀法向以及切向力的計算中ANALYN模型與KP模型均能得到一個較好的結果，而STRIPES模型在某些工況下會產生較大的誤差，不夠穩定，需要改進其曲率平滑方法以改善計算結果。相比于ANALYN模型，KP模型的計算過程簡單，計算時間短，因此只計算宏觀接觸力時建議采用KP模型。

為了能夠兼顧動力學計算與輪軌傷損研究的需要，這三種非赫茲接觸模型都需要進一步完善，以提高在不同接觸條件下的計算精度。需要指出的是，本文的分析僅針對標準的輪軌型面，而對于磨耗型面的影響則并未涉及，這部分工作有待在下一步研究中進行開展。