陳 雨，王攀杰，孫耀亮，安博洋，張 岷，王 平

(1.西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室，四川 成都 610031；2.西南交通大學(xué) 土木工程學(xué)院，四川 成都 610031；3.中鐵第一勘察設(shè)計(jì)院集團(tuán)有限公司，陜西 西安 710043)

由于鐵路運(yùn)營速度的提高和軸重的增加，輪軌接觸條件惡化，輪軌間磨耗加劇，這不僅會(huì)增加運(yùn)營成本，甚至還會(huì)危及行車安全。國內(nèi)外學(xué)者針對(duì)輪軌磨耗進(jìn)行了大量理論與試驗(yàn)研究，輪軌滾動(dòng)接觸模型作為其中重要的一環(huán)被應(yīng)用于輪軌磨耗仿真分析中。楊光等[1]采用Hertz接觸理論結(jié)合Fastsim算法研究了高速列車曲線通過時(shí)的輪軌磨耗情況；王璞等[2]研究了重載線路不同地段鋼軌磨耗的發(fā)展規(guī)律；安博洋等[3]采用基于虛擬滲透的非Hertz接觸簡(jiǎn)化算法分別結(jié)合Fastsim和FaStrip預(yù)測(cè)了接觸斑內(nèi)輪軌磨耗深度的分布；李霞等[4]采用Kalker三維非赫茲接觸模型預(yù)測(cè)了朔黃重載線路曲線部分的鋼軌磨耗，并與實(shí)測(cè)結(jié)果對(duì)比發(fā)現(xiàn)基本吻合，但是高軌的軌距角部分磨耗相差較大；孫宇等[5]采用Kalker三維非赫茲接觸模型預(yù)測(cè)鋼軌沿縱向和橫向磨耗的演變形態(tài)。以上均是對(duì)區(qū)間線路的分析，道岔區(qū)的接觸更加復(fù)雜。徐井芒等[6]采用半Hertz接觸算法結(jié)合FaStrip分析了軌道參數(shù)對(duì)道岔區(qū)鋼軌磨耗的影響；馬曉川等[7]對(duì)比分析了不同非赫茲滾動(dòng)接觸算法求解道岔區(qū)滾動(dòng)接觸行為時(shí)的計(jì)算精度和效率，建議道岔區(qū)采用Schani模型。

目前用于輪軌接觸分析的模型大都基于平面假設(shè)，但當(dāng)車輛經(jīng)過道岔或者小半徑曲線地段時(shí)，磨耗的輪緣與軌距角易發(fā)生共形接觸，此時(shí)的接觸斑是曲面形式，因此需要考慮更精確的模型計(jì)算輪軌接觸力，進(jìn)而分析輪軌磨耗和疲勞損傷等問題。

基于計(jì)算機(jī)和有限元模型的發(fā)展，不少學(xué)者利用有限元方法計(jì)算輪緣與軌距角發(fā)生共形接觸時(shí)的接觸應(yīng)力。金學(xué)松等[8]基于Kalker的三維彈性體非赫茲滾動(dòng)接觸理論，考慮曲面接觸斑，利用有限元法推導(dǎo)力與位移之間的關(guān)系，建立了任意曲面接觸斑狀態(tài)下的接觸余能表達(dá)式；李自力[9]利用有限元法計(jì)算準(zhǔn)四分之一空間的影響因子，替代Kalker三維滾動(dòng)接觸理論中的基于平面半空間假設(shè)的影響因子，可求解輪緣與軌距角的共形接觸問題；竺春海等[10]采用有限元參數(shù)二次規(guī)劃法求解輪軌共形接觸問題，得到了輪緣貼靠鋼軌軌距角的過程中接觸力的變化規(guī)律；Edwin等[11]基于李自力的研究，利用有限元得到輪緣與軌距角發(fā)生共形接觸時(shí)的影響因子，并與CONTACT數(shù)值程序結(jié)合，研究了共形接觸程度對(duì)接觸應(yīng)力分布的影響；楊新文等[12]利用有限元法計(jì)算輪軌法向接觸影響系數(shù)，對(duì)Kalker的三維非赫茲滾動(dòng)接觸理論中的法向應(yīng)力計(jì)算部分進(jìn)行修正，研究了重載鐵路曲線超高對(duì)鋼軌磨耗的影響。以上研究均利用有限元法得到輪軌間的影響因子進(jìn)而分析共形接觸問題，雖然提高了精度，但計(jì)算效率低，不適合大規(guī)模的計(jì)算。Blanco等[13]在Bossinesq-Cerruti半空間影響因子的基礎(chǔ)上進(jìn)行了修改，提出了一種求解任意空間影響因子的近似方法，并與有限元方法得到的影響因子進(jìn)行對(duì)比，結(jié)果基本吻合且效率得到了極大的提高；許玉德等[14]參考Blanco提出的影響因子，提出了最小余能方程中影響系數(shù)的修正公式，用于求解非平面接觸應(yīng)力。

現(xiàn)有文獻(xiàn)中未詳細(xì)對(duì)比分析基于平面假設(shè)的CONTACT數(shù)值程序與曲面算法的差異性，未指出在輪緣與軌距角發(fā)生共形接觸時(shí)，采用曲面算法的必要性。本文基于Kalker的三維非赫茲彈性體滾動(dòng)接觸理論(又稱“精確理論”)對(duì)其中的法向間隙、剛性蠕滑率和影響因子進(jìn)行修正，使其可考慮輪軌間曲面接觸斑，并對(duì)比分析了曲面算法和CONTACT在不同接觸情況下計(jì)算得到的接觸面積和接觸應(yīng)力等，研究了2種算法的一致性和差異性。本文的計(jì)算分析可為輪軌磨耗和疲勞損傷研究選擇適當(dāng)?shù)妮嗆墲L動(dòng)接觸模型提供參考。

1 考慮曲面接觸斑的輪軌滾動(dòng)接觸算法

車輛運(yùn)行過程中，會(huì)遇到小半徑曲線，這種地段輪緣與軌距角磨損通常較為嚴(yán)重，可能發(fā)生“共形接觸”。這種情況下，基于平面假設(shè)的Kalker精確理論只能得到近似解，基于此，本文對(duì)CONTACT的輸入?yún)?shù)進(jìn)行修改，使其適合求解共形接觸這種曲面形式的接觸斑問題。

1.1 基于平面假設(shè)的Kalker精確理論

Kalker精確理論通過將輪軌的潛在接觸區(qū)域離散化，基于彈性半空間假設(shè)，利用Boussinesq-Cerruti公式，得到三維彈性體非赫茲滾動(dòng)接觸最小余能的表達(dá)式[15-16]為

(1)

式中：AC為接觸斑區(qū)域；h為兩接觸體未變形前的法向間隙；δp為滲透量；un和uτ為接觸體在法向和切向的彈性位移差；pn和pτ為相應(yīng)方向的接觸應(yīng)力；u′τ表示前一時(shí)刻的切向彈性位移差；g=μpn為切向接觸邊界，μ為恒定摩擦系數(shù)；Wτ為切向剛性滑動(dòng)量；dx為接觸斑單元的縱向?qū)挾龋籨y為接觸斑單元的橫向?qū)挾取?/p>

彈性變形差u可以通過Boussinesq-Cerruti的力-位移方程求解

(2)

式中：上標(biāo)(1)、(2)分別為車輪和鋼軌；AIiJj為影響因子，表示當(dāng)單元J承受沿j方向的單位荷載時(shí)單元I處產(chǎn)生的沿i方向的彈性位移；AInJn為單元J承受沿法向的單位力時(shí)單元I處產(chǎn)生的沿法向的彈性位移；AIτJτ單元J承受沿切向的單位力時(shí)單元I產(chǎn)生的沿切向的彈性位移；Bij為j方向單位力引起車輪和鋼軌沿i方向的彈性位移差，其具體表達(dá)式可參考文獻(xiàn)[15]。

由于Kalker精確理論基于平面假設(shè)，適用于求解平面形式的接觸斑[16]，本文為了敘述方便，將其簡(jiǎn)稱為“平面算法”。

1.2 曲面算法

曲面算法與平面算法的主要區(qū)別在于考慮其接觸斑為曲面，接觸斑內(nèi)沿橫向分布的各點(diǎn)與接觸點(diǎn)的法向不再相同，這會(huì)影響接觸斑內(nèi)未變形前的法向間隙、剛性蠕滑率和影響因子的計(jì)算，而這些是計(jì)算輪軌法向和切向應(yīng)力的基本輸入?yún)⒘?，因此本文主要?duì)CONTACT的這3個(gè)參量進(jìn)行修改。

1.2.1 法向間隙

輪軌潛在接觸區(qū)域的法向間隙是計(jì)算輪軌法向接觸問題的基本幾何輸入變量。首先定義全局坐標(biāo)系OXYZ，原點(diǎn)O位于無位移時(shí)的輪對(duì)質(zhì)心處，X軸指向輪對(duì)前進(jìn)方向，Z軸垂直于軌道平面向上，Y軸水平向右，全局坐標(biāo)系固定不變。平面算法的接觸點(diǎn)坐標(biāo)系為Ocxyz，原點(diǎn)Oc位于輪軌接觸點(diǎn)，x軸指向車輪前進(jìn)方向，y軸指向Oc處的切線方向，z軸指向Oc處法線方向，接觸斑內(nèi)沿橫向y軸各點(diǎn)法向與接觸點(diǎn)的法向z一致；曲面算法的接觸點(diǎn)坐標(biāo)系為Ocxsn，x軸指向車輪前進(jìn)方向，s軸指向曲面接觸斑橫向，n軸指向接觸斑各點(diǎn)局部坐標(biāo)系法向，接觸斑沿橫向各點(diǎn)法向隨接觸面的幾何發(fā)生變化，見圖1。

圖1 2種算法的法向間隙示意圖

輪軌接觸斑內(nèi)質(zhì)點(diǎn)對(duì)產(chǎn)生接觸的條件是接觸質(zhì)點(diǎn)對(duì)到接觸點(diǎn)的弧長(zhǎng)相等。以接觸斑內(nèi)E點(diǎn)為例說明2種算法計(jì)算法向間隙的區(qū)別，見圖1，BE與接觸點(diǎn)處法向Ocz平行，平面算法中考慮BE段的距離即為輪軌間法向間隙，這在接觸點(diǎn)附近輪軌曲率半徑很大時(shí)可以滿足；但是在輪緣與軌距角形成曲面接觸的情況下，BOc與EOc距離并不相等，即此時(shí)B、E并未產(chǎn)生接觸，因此平面算法會(huì)高估真實(shí)的法向間隙。

曲面算法求解法向間隙的思想在文獻(xiàn)[11，13]中有提及，但是未給出數(shù)學(xué)表達(dá)式，本文采用一種近似方法。根據(jù)切平面坐標(biāo)系Ocxyz和曲面坐標(biāo)系Ocxsn可建立起(y，z，s)之間的對(duì)應(yīng)關(guān)系，OcE和OcC弧長(zhǎng)坐標(biāo)sw和sr相等，那么接觸斑內(nèi)D點(diǎn)對(duì)應(yīng)的法向間隙ds為

(3)

式中：θ為EC間距離與D點(diǎn)的局部法向n間的夾角。

分別采用平面算法和曲面算法得到輪軌二維法向間隙f(y)后，假設(shè)沿縱向x方向仍遵循Hertz橢圓分布，可得三維的法向間隙

h=Ax2+f(y)

(4)

式中：x為接觸斑內(nèi)縱坐標(biāo)；A=1/(2R)，平面接觸算法中R=R0恒為接觸點(diǎn)的滾動(dòng)圓半徑，曲面接觸算法中R是接觸曲面沿s軸各點(diǎn)的滾動(dòng)圓半徑。

1.2.2 剛性蠕滑率

輪軌潛在接觸區(qū)域內(nèi)的剛性蠕滑率是產(chǎn)生輪軌切應(yīng)力的主要因素。根據(jù)蠕滑率定義，剛性蠕滑率為

式中：wY為車輪繞Y軸旋轉(zhuǎn)的角速度；δ為接觸斑各點(diǎn)的接觸角，見圖2；v0為輪軌接觸界面的平均速度，鋼軌靜止不動(dòng)，可假設(shè)v0=wYR0；εx,εy分別為接觸點(diǎn)處的縱、橫向蠕滑率。

圖2 局部坐標(biāo)系

參考文獻(xiàn)[9]，基于平面假設(shè)，接觸斑內(nèi)各點(diǎn)接觸角恒等于接觸點(diǎn)處接觸角，則R0-R=ysinδ0，xsinδ=xsinδ0，平面算法的剛性滑移率為

(7)

曲面算法中，由于接觸斑是曲面，接觸角并不恒等于δ0，因此各點(diǎn)的剛性蠕滑率為

(8)

在曲面接觸算法中，取接觸斑內(nèi)沿橫向各點(diǎn)的接觸角δ=δ0時(shí)，則可得到平面算法剛性蠕滑率公式。

1.2.3 彈性位移差

輪軌潛在接觸區(qū)域內(nèi)的影響因子是求解法向和切向接觸問題的基本輸入?yún)⒘?。在以往的研究[9,11]中，通常采用有限元方法獲得任意曲面空間下的彈性位移差，但這會(huì)顯著地增加計(jì)算成本，難以應(yīng)用于工程計(jì)算。為了提高計(jì)算效率，本文采用Blanco等[13]提出的可考慮任意曲面的彈性位移差計(jì)算方法(該方法的有效性已通過與有限元方法的對(duì)比得到了驗(yàn)證)。曲面接觸算法求解彈性位移差的示意圖圖3，α為接觸斑內(nèi)I點(diǎn)與J點(diǎn)的切向角之差。在平面算法中，接觸斑內(nèi)各點(diǎn)的法向與接觸點(diǎn)的法向一致，α=0；在曲面算法中，考慮整個(gè)接觸斑內(nèi)法向沿接觸曲面不斷變化，α≠0。

圖3 曲面算法彈性位移差的示意

曲面接觸的彈性位移差為

(9)

式中：Bij、Bij_c分別為平面算法、曲面算法的彈性位移差，其中,i,j分別為切平面坐標(biāo)系中的x、y、z方向、曲面坐標(biāo)系中的x、s、n方向。在平面算法中，對(duì)于同種材料，法向與切向互不影響，即Byn、Bxn、Bnx、Bny均為0；在曲面算法中，考慮了法向與切向的耦合。

2 考慮曲面接觸斑的輪軌接觸行為分析

為了對(duì)比平面算法與曲面算法研究輪軌接觸行為的不同，本文以中國S1002CN車輪踏面和CHN60鋼軌型面為例進(jìn)行計(jì)算。車輪名義滾動(dòng)圓半徑為460 mm，軌距為1 435 mm，軌底坡為1/40。輪軌接觸特性的計(jì)算參數(shù)：輪軌垂向力為83.3 kN，輪軌摩擦系數(shù)為0.3，泊松比為0.28，剪切模量為82 GPa。輪軌間接觸位置大致可以分為3種：車輪踏面與軌頂接觸、車輪輪緣與軌距角共形接觸和車輪輪緣與軌距角凸面接觸，本文主要選取了3個(gè)不同的輪對(duì)橫移量，分別為-6、7.77、9.5 mm代表這3種接觸位置。下面主要對(duì)這3種典型的輪軌接觸情況進(jìn)行詳細(xì)分析。

2.1 輪軌法向接觸行為分析

分別采用平面算法和曲面算法研究輪軌法向接觸行為，計(jì)算內(nèi)容包括輪軌法向間隙、接觸斑面積以及法向接觸應(yīng)力分布。

2.1.1 輪軌法向間隙分布

輪軌間的法向間隙會(huì)直接影響到接觸斑面積和接觸應(yīng)力，分別采用平面算法和曲面算法計(jì)算3種接觸情況下的輪軌法向間隙，得到的結(jié)果見圖4，左圖為輪軌接觸點(diǎn)在全局坐標(biāo)系下的位置，右圖為輪軌法向間隙在接觸點(diǎn)坐標(biāo)系下沿橫向的分布。由圖4可見，當(dāng)輪對(duì)橫移量為-6 mm時(shí)，車輪踏面與軌頂接觸，輪軌型面在接觸點(diǎn)附近的曲率半徑與接觸斑幾何尺寸相比很大，2種算法得到的法向間隙幾乎沒有差異；當(dāng)輪軌橫移量為7.77 mm時(shí)，輪緣與軌距角發(fā)生共形接觸，輪軌在接觸點(diǎn)附近曲率半徑均較小且曲率中心在同側(cè)，這種情況下平面算法得到的法向間隙明顯比曲面算法大。這是由于曲面算法考慮輪軌接觸斑的真實(shí)法向，而平面算法考慮接觸斑內(nèi)法向均與接觸點(diǎn)處的法向相同；在輪緣與軌距角發(fā)生共形接觸時(shí)，接觸斑內(nèi)法向差異較大，因此平面算法會(huì)高估其真實(shí)的法向間隙；當(dāng)輪對(duì)橫移量為9.5 mm時(shí)，輪緣與軌距角發(fā)生凸面接觸，這種情況雖然輪軌型面在接觸點(diǎn)附近曲率半徑小且接近，但是曲率中心不在同一側(cè)，可以將輪軌接觸面近似看成平面，2種算法得到的法向間隙也幾乎沒有差異。

圖4 不同橫移量下的輪軌接觸點(diǎn)位置和接觸點(diǎn)坐標(biāo)系中法向間隙沿橫向分布

2.1.2 輪軌接觸斑面積和法向應(yīng)力分布

由于輪對(duì)橫移量為7.77 mm時(shí)2種算法在接觸點(diǎn)坐標(biāo)系內(nèi)得到的法向間隙差別很大，因此本節(jié)分別采用平面算法和曲面算法進(jìn)一步計(jì)算這種接觸情況下的接觸斑和法向應(yīng)力分布，結(jié)果見圖5。由圖5(a)可知，曲面算法得到的左接觸斑面積明顯大于平面算法，其接觸斑面積約為84 mm2，平面算法約為76 mm2，比曲面算法的結(jié)果小9.5%。原因是平面算法得到的法向間隙更大，因此在相同的法向力作用下其接觸面積更小。由圖5(b)可見，左接觸斑內(nèi)，曲面算法得到的最大法向應(yīng)力約為1 732 MPa，平面算法約為2 328 MPa，比曲面算法的結(jié)果大34.4%。說明當(dāng)輪緣與軌距角發(fā)生共形接觸時(shí)，平面算法會(huì)高估真實(shí)情況下的法向應(yīng)力。

圖5 橫移量為7.77 mm時(shí)在接觸點(diǎn)坐標(biāo)系下的接觸斑和法向應(yīng)力

2.2 輪軌切向接觸行為分析

研究輪軌切向接觸行為的主要目的是確定輪軌切應(yīng)力分布，以下將分別采用平面算法和曲面算法分析輪軌切向接觸行為。

2.2.1 輪軌接觸斑內(nèi)接觸角變化

輪軌接觸斑內(nèi)的接觸角會(huì)影響自旋，從而導(dǎo)致縱橫向剛性蠕滑率變化。平面算法考慮接觸斑內(nèi)的接觸角恒等于接觸點(diǎn)處的接觸角，但實(shí)際上，接觸斑內(nèi)接觸角會(huì)隨著輪軌型面發(fā)生變化，尤其是輪緣與軌距角發(fā)生曲面接觸時(shí)。分析輪對(duì)橫移量分別為-6、7.77、9.5 mm 3種接觸情況的接觸角在接觸斑內(nèi)的變化情況，見圖6。

圖6 不同橫移量下的接觸角

由圖6可知，輪對(duì)橫移量為-6 mm時(shí)，接觸斑寬度約為11.2 mm，接觸斑內(nèi)接觸角變化幅度約為1.24°；輪對(duì)橫移量為9.5 mm時(shí)，接觸斑寬度約為0.85 mm，此時(shí)接觸斑沿縱軸呈狹長(zhǎng)狀，接觸斑內(nèi)接觸角變化幅度約為0.04°。踏面與軌頂發(fā)生接觸或者輪緣與軌距角呈凸面接觸時(shí)，變化幅度均較小，此時(shí)平面算法將接觸斑內(nèi)的接觸角假設(shè)為與接觸點(diǎn)處的接觸角一致，不會(huì)產(chǎn)生較大誤差。而輪對(duì)橫移量為7.77 mm時(shí),接觸斑寬度約為11.1 mm，接觸斑內(nèi)接觸角變化幅度達(dá)到了22.68°。根據(jù)圖2和式(7)、式(8)知，接觸角變化會(huì)影響接觸斑內(nèi)自旋，曲面算法考慮了真實(shí)接觸角的變化，在輪緣與軌距角發(fā)生共形接觸時(shí)，其縱橫向蠕滑率會(huì)與平面算法相差較大，進(jìn)而影響到切向力的大小和分布。

2.2.2 輪軌影響因子

圖7 不同曲率下鋼軌和車輪的影響因子

由圖7可知，車輪和鋼軌曲率半徑為300 mm時(shí)，曲面算法得到的車輪和鋼軌變形量Asn幾乎與平面算法完全一致；半徑為80 mm時(shí)有輕微差別；然而半徑為13 mm時(shí)差別較明顯。原因是半徑為300 mm時(shí)，接觸點(diǎn)附近接觸型面法向差異不大，即式(9)中α角小，導(dǎo)致曲面算法得到的影響因子與平面算法差異較小；當(dāng)半徑為13 mm時(shí)，α角大，從而導(dǎo)致2種算法得到的影響因子差異較大。此外，在相同的法向力下，車輪和鋼軌的橫向位移并不相同，輪軌間的相對(duì)位移會(huì)產(chǎn)生蠕滑現(xiàn)象，從而產(chǎn)生切應(yīng)力。

2.2.3 蠕滑曲線

車輛動(dòng)力學(xué)中，輪軌力是重要的參量，求解蠕滑力的基礎(chǔ)是蠕滑率、蠕滑力關(guān)系曲線。由于3個(gè)工況中，輪對(duì)橫移量為7.77 mm時(shí)接觸斑內(nèi)接觸角變化最大，因此本節(jié)選擇輪軌橫移量為7.77 mm時(shí)的蠕滑率、蠕滑力關(guān)系曲線進(jìn)行分析。為了得到曲面算法由于考慮變自旋所導(dǎo)致的蠕滑率、蠕滑力關(guān)系曲線與平面算法的差異，分別采用平面算法和曲面算法計(jì)算輪軌接觸斑內(nèi)的蠕滑率、蠕滑力關(guān)系曲線，僅改變接觸點(diǎn)處的縱向或橫向蠕滑率，即式(7)、式(8)中的εx或εy，縱坐標(biāo)除以切向力邊界10-6P使其無量綱化，結(jié)果見圖8。圖8(a)為縱向蠕滑率改變而橫向蠕滑率為0的蠕滑曲線，圖8(b)為橫向蠕滑率改變而縱向蠕滑率為0的蠕滑曲線。

圖8 蠕滑率、蠕滑力關(guān)系曲線對(duì)比

由圖8可知，2種算法得到的蠕滑率和蠕滑力關(guān)系均表現(xiàn)出明顯的非線性特征，這是由于蠕滑力受到庫倫摩擦定律限制，當(dāng)蠕滑率達(dá)到一定值時(shí)，接觸斑進(jìn)入全滑動(dòng)狀態(tài)，蠕滑力達(dá)到極限值。圖8平面算法得到的縱向蠕滑力和橫向蠕滑力在橫坐標(biāo)為0時(shí)縱坐標(biāo)不為0是由于考慮了自旋的影響，且相對(duì)誤差均在橫坐標(biāo)為0時(shí)最大，分別為9.1%和9.4%。2種算法得到的蠕滑曲線不同的原因是曲面算法考慮了接觸斑內(nèi)變自旋，即接觸斑內(nèi)接觸角和滾動(dòng)圓半徑變化，而平面算法考慮恒定自旋，從而引起縱橫向剛性蠕滑率的差異，蠕滑率、蠕滑力曲線的不同則會(huì)導(dǎo)致蠕滑力的不同。

2.2.4 輪軌切向接觸應(yīng)力大小及分布

分別采用平面算法和曲面算法計(jì)算橫移量為-6、7.77、9.5 mm時(shí)的切向力，結(jié)果見圖9。圖9中，Xc為接觸斑內(nèi)沿縱向的坐標(biāo)，Yc為接觸斑內(nèi)沿橫向的坐標(biāo)，顏色深淺代表切應(yīng)力的大小，箭頭代表切應(yīng)力的方向。由圖9(a)、9(b)、9(e)、9(f)可知，當(dāng)輪對(duì)橫移量為-6、9.5 mm時(shí)，平面算法和與曲面算法得到的切應(yīng)力的大小和方向幾乎沒有差異。原因是這2種接觸情況下，接觸斑內(nèi)接觸角變化小，2種算法得到的縱橫向蠕滑率差異較小。雖然橫移量為9.5 mm時(shí)，輪緣與軌距角接觸，輪軌型面接觸半徑大約為13 mm，曲面算法得到的影響因子與平面算法得到的結(jié)果差別較大，但是這種差別對(duì)接觸應(yīng)力的影響不大。由圖9(c)、9(d)可知，當(dāng)輪軌橫移量為7.77 mm時(shí)，平面算法和曲面算法得到的切應(yīng)力大小明顯不同。曲面算法得到的最大切應(yīng)力約為698 MPa，平面算法約為519 MPa，比曲面算法的結(jié)果大34.5%。原因是橫移量為7.77 mm時(shí)，輪緣與軌距角發(fā)生共形接觸，這時(shí)2種算法得到的影響因子差別較大,且此時(shí)接觸斑內(nèi)接觸角變化大，從而使2種算法得到的縱橫向蠕滑率差別也較大?；谝陨系姆治隹芍绊懸蜃雍腿浠实木C合影響導(dǎo)致2種算法的切應(yīng)力結(jié)果差異較大，且主要是由于蠕滑率的影響。

圖9 不同橫移量下接觸斑內(nèi)切應(yīng)力分布

為了更直觀地體現(xiàn)出輪緣與軌距角發(fā)生曲面接觸時(shí)2種算法得到的切應(yīng)力分布差別，分析接觸斑橫坐標(biāo)Yc為-7 mm斷面上接觸切應(yīng)力沿接觸斑縱向分布，結(jié)果見圖10。由圖可知2種算法得到的縱橫向切應(yīng)力差異較大。曲面、平面算法得到的最大縱橫向切應(yīng)力分別為257、36 MPa；平面算法分別約為452、78 MPa，平面算法得到的最大縱橫向切應(yīng)力分別比曲面算法的結(jié)果大75.9%、113.5%。由此可以看出2種算法的切應(yīng)力大小及分布差異較大，這種差異會(huì)對(duì)輪軌磨耗深度的計(jì)算產(chǎn)生較大影響。

圖10 輪對(duì)橫移量為7.77 mm時(shí)的縱、橫向切應(yīng)力

2.3 輪軌接觸斑內(nèi)磨耗深度分布

分析發(fā)現(xiàn)2種算法在輪緣與軌距角發(fā)生共形接觸時(shí)計(jì)算得到的輪軌接觸行為差異較為明顯，且車輪過小半徑曲線時(shí)，輪緣與軌距角易發(fā)生共形接觸從而加劇輪軌磨耗，因此本節(jié)采用2種算法分析輪緣與軌距角發(fā)生共形接觸時(shí)的磨耗深度。利用平面算法和曲面算法計(jì)算得到輪軌間蠕滑率和切應(yīng)力，然后采用USFD磨耗模型[17]預(yù)測(cè)磨耗深度，由圖5分析可知，輪對(duì)橫移量為7.77 mm時(shí)輪緣與軌距角會(huì)發(fā)生共形接觸，因此本節(jié)以橫移量為7.77 mm為例分析2種算法預(yù)測(cè)磨耗深度的差異。為了便于對(duì)比，將2種算法得到的磨耗深度除以平面算法的最大磨耗深度，進(jìn)行歸一化處理，沿接觸斑橫向坐標(biāo)磨耗深度的結(jié)果見圖11。

圖11 橫移量為7.77 mm時(shí)的預(yù)測(cè)磨耗深度

由圖11可知，2種算法得到的最大磨耗深度不同，曲面算法得到的最大磨耗深度約為平面算法的93.5%；此外，2種算法得到的左側(cè)接觸斑磨耗范圍也不同，平面算法的磨耗范圍為-8.6～-5.9 mm，曲面算法磨耗范圍為-9.5～-6 mm。其原因主要是2種算法得到的接觸斑不同、切應(yīng)力和蠕滑率大小和分布不同?？梢姰?dāng)輪緣與軌距角發(fā)生共形接觸時(shí)，2種算法得到的磨耗深度差異較明顯。因此，建議采用曲面算法計(jì)算磨耗深度，以便更加精確地預(yù)測(cè)鋼軌磨耗的演變。

2.4 考慮輪軌曲面接觸行為的必要性分析

通過分析可知，輪緣與軌距角發(fā)生共形接觸時(shí)平面算法的結(jié)果會(huì)產(chǎn)生較大誤差，因此有必要采用曲面算法。雖然上文的分析中只選取了輪對(duì)橫移量為7.77 mm這一個(gè)共形接觸工況，但這不僅僅是特例。為了說明這個(gè)問題，分別采用平面算法和曲面算法計(jì)算了輪對(duì)不同橫移量下的接觸斑面積和最大法向應(yīng)力，橫移量從-6～12 mm，步進(jìn)為0.5 mm，在7.77 mm附近加密取了7.76、7.78、7.79 mm，結(jié)果見圖12。

圖12 不同橫移量下的接觸斑面積和最大法向應(yīng)力及平面算法的相對(duì)誤差

由圖12(a)、12(c)可知，隨著橫移量的增加，接觸斑面積大致呈減小趨勢(shì)，最大法向應(yīng)力大致呈增大趨勢(shì)。但是在7.77 mm附近，接觸斑面積先增大后減小，最大法向應(yīng)力先減小后增大，這是由于輪緣與軌距角發(fā)生共形接觸時(shí)接觸斑面積會(huì)增大。此外，由圖12(b)、12(d)可知，平面算法相對(duì)誤差均集中在7.77 mm附近分布，接觸斑面積相對(duì)誤差最大約為13%，出現(xiàn)在橫移量為7.79 mm處；最大法向應(yīng)力相對(duì)誤差最大為35.6%，出現(xiàn)在橫移量為7.78 mm處。

綜上所述，輪緣與軌距角發(fā)生共形接觸時(shí)，接觸斑形式是曲面，此時(shí)采用平面算法求解輪軌接觸會(huì)產(chǎn)生較大的誤差，為了更精確地分析輪緣與軌距角發(fā)生共形接觸時(shí)的輪軌滾動(dòng)接觸行為，有必要采用曲面算法。

3 結(jié)論

為了更精確地計(jì)算輪緣與軌距角發(fā)生共形接觸時(shí)的接觸斑面積、接觸應(yīng)力，本文基于Kalker三維非赫茲滾動(dòng)接觸理論，對(duì)CONTACT數(shù)值程序的法向間隙、縱橫向蠕滑率和影響因子進(jìn)行修改使其可以考慮輪軌曲面接觸斑。以S1002CN車輪踏面和CHN60鋼軌型面為研究對(duì)象，對(duì)比分析了曲面算法與CONTACT在踏面與軌頂發(fā)生接觸、輪緣與軌距角發(fā)生共形接觸、輪緣與軌距角發(fā)生凸面接觸這3種典型接觸工況下的計(jì)算結(jié)果，得到如下結(jié)論。

(1)當(dāng)車輪踏面與軌頂接觸或者輪緣與軌距角呈凸面接觸時(shí)，接觸斑內(nèi)沿橫向各點(diǎn)法向差異較小，可以將接觸斑看成平面，此時(shí)平面算法和曲面算法得到的法向間隙、接觸斑面積和接觸應(yīng)力差異均較小。

(2)當(dāng)車輪輪緣與軌距角發(fā)生共形接觸時(shí)，接觸斑內(nèi)沿橫向各點(diǎn)法向差異較大，導(dǎo)致采用2種算法得到的輪軌滾動(dòng)接觸行為差異也較大。主要體現(xiàn)在：平面算法得到的法向間隙比曲面算法得到的結(jié)果大，接觸斑面積比曲面算法得到的結(jié)果小。因此，平面算法的最大法向接觸應(yīng)力和切向接觸應(yīng)力均比曲面算法得到的結(jié)果大。此外，2種算法得到的蠕滑率、蠕滑力關(guān)系曲線也不同。在該接觸情況下，使用平面算法得到的結(jié)果會(huì)產(chǎn)生較大的誤差，因此有必要采用曲面算法。

(3)將2種算法得到的蠕滑率和切應(yīng)力用于計(jì)算輪緣與軌距角發(fā)生共形接觸時(shí)的磨耗分布，發(fā)現(xiàn)其計(jì)算結(jié)果差異較大。因此建議這種接觸情況下采用曲面算法計(jì)算磨耗，以便更精確地預(yù)測(cè)輪軌磨耗演變。

本文通過理論推導(dǎo)，提出了一種基于CONTACT邊界元方法的曲面接觸模型，初步調(diào)查了曲面接觸特性對(duì)滾動(dòng)接觸解的影響，后續(xù)將進(jìn)一步通過有限元方法對(duì)其有效性進(jìn)行驗(yàn)證，并將其植入車輛-軌道耦合動(dòng)力學(xué)程序，考慮曲面接觸特性對(duì)車輛動(dòng)力學(xué)行為的影響。