李糧余 田春香, 周佳儀 安博洋 王 平

(1.中鐵二院工程集團有限責任公司,成都 610031；2.西南交通大學高速鐵路線路工程教育部重點實驗室,成都 610031)

1 研究背景

高速鐵路輪軌關系是制約列車運行速度以及影響車輛安全運營的關鍵問題[1-2]，因此，需要有效的輪軌接觸模型以分析輪軌間的相互作用，能夠準確、快速地求解非赫茲滾動接觸行為。由于基于邊界元法[3]和有限元法[4]的準確接觸模型的計算時間較長，輪軌磨耗預測需要進行數(shù)百萬次計算以模擬實際的運行距離，因此，無法高效預測輪軌磨耗。

在提高計算效率的同時考慮輪軌非赫茲滾動接觸行為，一些學者提出了幾種簡化的非赫茲接觸模型。Kik和Piotrowski[5]提出了“虛擬滲透”的概念，假設壓應力沿橫向的變化與接觸斑縱向邊界的分布保持一致，該模型可快速地確定接觸斑形狀?；谔摂M滲透方，Linder提出了另一種接觸斑形狀的確定方法，即將接觸斑視為一系列沿著縱向分布的條帶。然而，Kik和Linder等人的方法均無法獲得較為準確的壓應力分布。相比較而言，Ayasse和Chollet[6]開發(fā)的STRIPES接觸模型，采用赫茲接觸理論求解每一條帶上的壓應力分布，可獲得較為準確的結果。類似地，Sichani[7]等提出了一種名為ANALYN的接觸模型，考慮了彈性變形對接觸斑和應力的影響。

上述簡化接觸模型的切向計算一般采用Kalker簡化理論及其FASTSIM算法，由于FASTSIM是基于橢圓接觸假設而開發(fā)的，因此，處理非赫茲接觸問題時需對橢圓計算所得到的柔度系數(shù)進行修正。在Kik和Piotrowski的模型中，通過接觸斑面積和長短半軸比兩個參數(shù)定義一個等效橢圓，進而計算該橢圓的柔度系數(shù)，并作為非赫茲接觸斑的切向計算參數(shù)；Linder和STRIPES接觸模型則將接觸斑內每個條帶等效為一個局部橢圓，并計算相應的柔度系數(shù)。盡管按照上述方法可得到較為準確的蠕滑力，但由于FASTSIM假設接觸應力與彈性位移呈線性關系并且引入了不真實的拋物線型切向應力邊界條件，這些方法仍無法獲得較為準確的切向接觸結果。為此，Sichani提出了一種基于Kalker的條帶理論和簡化理論的新型接觸計算方法——FaStrip，可獲得更為準確的切向接觸結果。

由于STRIPES接觸模型可較好地兼顧計算精度與效率，已廣泛應用于車輛軌道/道岔動力學、接觸力學建模以及磨耗預測。但在最近的研究中發(fā)現(xiàn) STRIPES接觸模型在某些極端非赫茲工況下存在較大計算誤差[8]。因此，本文對STRIPES接觸模型進行改進，提高其計算精度和穩(wěn)定性。

2 STRIPES算法

2.1 法向接觸計算

輪軌法向接觸問題中，法向間隙f(y)、滲透量δ以及彈性變形u(y)之間存在以下關系：

(1)

式中：R——接觸點處車輪半徑(m);

x——滾動方向的坐標；

y——橫向方向的坐標。

應力p(x,y)是在物體表面產生的彈性變形，可用Boussinesq-Cerruti方程表示：

(2)

式中：Aij——影響系數(shù)，代表單元j處的單位載荷引起單元i處的變形；

V——材料的泊松比；

E——材料的彈性模量(MPa)。

為簡化求解式(1)的計算過程，STRIPES基于虛擬滲透原理，采用比例系數(shù)ε替代計算彈性變形。在該算法中，將接觸斑沿滾動方向劃分為條帶，如圖1所示。接觸斑的橫向邊界可簡化為：

h(y)=εδ-f(y)≥0

(3)

式中，ε由赫茲理論計算得到：

(4)

式中：A0——第一個接觸點處的縱向曲率；

B0——第一個接觸點處的橫向曲率；

n0、r0——無量綱赫茲常數(shù)。

假設每個條帶對應一個局部等效橢圓，且具有相同的赫茲特性，每個條帶的半軸長度a(y)可利用式(5)求解：

(5)

式中：m——無量綱赫茲常數(shù)。

接觸斑內任意單元的應力為：

(6)

2.2 切向接觸計算

在STRIPES中切向接觸問題采用Kalker的FASTSIM算法進行處理，假設每個條帶具有各自的橢圓特性，在每個條帶中使用基于橢圓接觸開發(fā)的FASTSIM算法進行計算，如圖1所示。在剪切應力設定為0的接觸斑前沿計算每個條帶的剪切應力增量：

圖1 STRIPES模型示意圖

(7)

(8)

式中：υx(y)——每個條帶處縱向蠕滑率；

φ(y)——每個條帶處的自旋；

υy——橫向蠕滑率；

L1、L2、L3——每個局部橢圓的柔度系數(shù)。

(9)

式中：c11,c22,c23——Kalker蠕滑力系數(shù)；

G——剪切模量(MPa)。

假設牽引邊界遵循拋物線分布：

(10)

式中：μ——摩擦力系數(shù)。

(11)

否則，該點位于滑動區(qū)，此時剪切應力為：

(12)

在實際情況中牽引力界限并非拋物線型，但這一假設對于確保在FASTSIM中得到相對準確的蠕滑力曲線是必不可少的，因此，在實際切向計算中，牽引力邊界由式(6)中橢圓分布變?yōu)槭?10)中的拋物線分布。有關此模型的更多詳細信息，請參閱文獻[7]。

3 STRIPES模型修正

3.1 由接觸邊界確定可變比例系數(shù)

從式(4)可以看出，STRIPES中假定比例系數(shù)ε為常量，忽略了接觸區(qū)域上曲率的變化，接觸斑形狀的確定僅取決于第一接觸點。因此，本文通過采用每個條帶處變化的比例系數(shù)ε(y)替換單一系數(shù)ε從而對STRIPES算法進行修正：

(13)

將式(13)代入式(3)、式(5)和式(6)，可獲得新的接觸斑區(qū)域以及應力分布。與KP和Linder等其他虛擬滲透方法相比，STRIPES接觸模型的優(yōu)點是當接觸斑為橢圓形的情況下可得到與赫茲理論一致的結果，而其他兩種方法則不能。由于橢圓接觸時曲率保持恒定，考慮可變比例系數(shù)并不會改變STRIPES接觸算法的特性，此時ε(y)保持恒定且等于ε，因此ε屬于ε(y)的一種特殊情況。

3.2 一種獲得合理滲透量的簡單迭代方法

與其他接觸模型一樣，STRIPES模型式(3)中的滲透量δ也是未知的，為使模型計算所得法向力與給定法向力一致，需迭代求解。之前的研究中常使用赫茲理論所得滲透量作為δ，原因在于LMa-CN60組合輪軌廓形產生的接觸斑接近于橢圓形。然而在非橢圓形的接觸情況下赫茲滲透量和實際滲透量之間可能存在較大差異，因此有必要增加迭代過程以獲得真實滲透量。據此，本節(jié)提出了一種實現(xiàn)該過程的新方法。

以給定的法向力N為已知參數(shù)，由式(14)計算所得赫茲滲透量δH作為初始滲透量。

(14)

按式(15)進行迭代，直至由式(6)應力積分計算得到的輸出法向力P與給定的法向力N相等。

δn=λδn-1

(15)

式中：λ——無量綱比例因子。

(16)

由于法向力對滲透量值的變化非常敏感，較大的比例因子可能使輸出法向力產生劇烈波動，而小值導致不能快速收斂。式(16)中所定數(shù)值是由大量工況計算確定以確保能夠加速迭代過程。盡管需要進行滲透量迭代，簡化模型的計算速度仍比精確模型快得多。

3.3 與FaStrip相結合的非橢圓方法

如2.2節(jié)所述，STRIPES借助Kalker的FASTSIM算法進行切向計算，然而該算法由于假設剪切應力在黏著區(qū)呈線性分布且遵循拋物線牽引力邊界而不能得到令人滿意的滾動接觸解。為提高計算精度，本節(jié)引入名為FaStrip的新切向接觸模型來代替FASTSIM。

FaStrip算法將Kalker條帶理論、線性理論與FASTSIM算法結合起來，可以考慮任意組合的縱向、橫向以及自旋蠕滑率。該方法假設接觸區(qū)為橢圓形并且沿滾動方向將其離散成條帶，通過改進的條帶理論計算黏著區(qū)的剪切應力分布。為保持與線性理論相同的蠕滑曲線斜率，該理論引入了3個校正因子，且滑動區(qū)的剪切應力分布遵循橢圓牽引界限。因此，F(xiàn)aStrip算法可在任意接觸橢圓形狀的情況下提供足夠準確的剪切應力和滑動量值，并保持與FASTSIM相同的計算效率。有關FaStrip的更多詳細信息，請參閱文獻[7]。

此外，為了適應非橢圓接觸條件，基于橢圓假設開發(fā)的FaStrip算法應對由橢圓條件計算的條帶柔度系數(shù)進行修正。這種改進可以確?；茀^(qū)域中剪切應力方向的正確性，否則將影響蠕滑力的橫縱向分量。

FaStrip算法中仍然使用式(9)中每個條帶的3個柔度系數(shù)，將式(7)修正為：

(17)

此改進充分考慮了自旋對縱向剪切應力的影響，當車輪與軌頂相接觸時，接觸角的變化較小，輪軌接觸遵循半空間假設，接觸斑可被視為平面，接觸斑內自旋和蠕滑也都保持不變。

另一方面，滑動區(qū)飽和剪切應力受限于式(18)中的橢圓牽引邊界：

(18)

4 計算結果分析

修正STRIPES的接觸模型(Modified STRIPES)實現(xiàn)了上述改進，該算法基于MATLAB軟件開發(fā)，計算速度約比FORTRAN程序編碼的Kalker 的CONTACT程序快200～300倍。

本節(jié)旨在評估修正STRIPES接觸模型計算法向和切向接觸的能力。采用軌底坡1∶40的輪軌標準型面組合S1002CN-CN60，輪軌界面之間的摩擦系數(shù)為0.3。-3～7 mm 10種不同橫向位移下接觸形狀和壓力的變化情況，如圖2所示。本文取輪對放置在中心位置(Δy= 0 mm)且橫向位移分別為2 mm和 6 mm時的工況作為典型接觸工況進行詳細分析。

圖2 CONTACT預測的橫移量范圍為-3～6 mm(從左到右)時接觸形狀、壓力變化圖

4.1 法向接觸

以Kalker的CONTACT程序作為對照，修正STRIPES接觸模型與原始STRIPES模型計算的沿橫向(在x=0處)接觸斑和法向應力分布情況進行了比較，如圖3～圖5所示。從圖3～圖5中可以看出，Δy=0 mm時，接觸斑處于高度非赫茲情況，右邊界大于左邊長度。與STRIPES相比，修正STRIPES模型計算的接觸斑和最大應力結果均得到一定的改善，接觸斑寬度和最大應力的誤差分別從22.06%降至8.94%、12.57%降至4.13%。Δy=2 mm時，三種接觸模型中的接觸斑和應力分布都比較接近，相較之下修正STRIPES模型計算結果更為準確，與CONTACT結果基本一致。Δy=6 mm,時接觸斑由兩個橢圓合并組成，STRIPES無法正確估計右側的真實接觸形狀，而修正后的STRIPES模型明顯改善了接觸斑形狀和最大應力的結果，更接近參考值。

圖3 三種工況接觸斑圖(x=0)

圖4 三種工況應力分布圖(x=0)

圖5 三種工況幾何接觸間隙圖(x=0)

上述改進歸因于修正STRIPES模型中使用的變化比例系數(shù)的影響。通過對不同橫移量時幾何間隙以及修正STRIPES模型與原始STRIPES模型中分別使用的滲透量值ε(y)δ、εδ對比發(fā)現(xiàn)，由第一接觸點估計的εδ能夠很好地計算左側邊界，原因在于幾何間隙的左邊部分接近橢圓接觸間隙。因此，三種接觸模型在接觸斑形狀和應力的負方向上能夠相互對應。而此時幾何間隙的右側顯示非橢圓特征，STRIPES模型會低估接觸區(qū)域的正橫向邊界。相比之下，使用變化比例系數(shù)可明顯改善法向接觸解。值得一提的是，在Δy= 6 mm工況下，兩種簡化模型比例系數(shù)與滲透量的乘積并不相同，滲透量會因第3.3節(jié)中所述的迭代過程而發(fā)生變化。

4.2 切向接觸結果

三種典型工況下的切向接觸結果如圖6所示，箭頭指向剪切應力的方向，長度與剪切應力的大小成正比。從計算結果中可發(fā)現(xiàn)，兩種簡化模型都可獲得令人比較滿意的黏滑分布和剪切應力值，其中修正STRIPES模型結果與CONTACT更為接近。

圖6 不同橫移量下切向力分布對比圖

但兩種簡化接觸模型預測的剪切應力分布在某些工況下仍存在一定差異。前兩種工況下，接觸區(qū)域中黏著區(qū)域居多，STRIPES模型預測的剪切應力明顯小于其他兩種算法，而最后一種工況則相反。這種現(xiàn)象主要是由于采用不同的牽引力邊界條件以及STRIPES模型在黏著區(qū)域中假設剪切應力呈線性分布所導致。三種典型工況的總剪切應力及其分量沿縱向的分布情況(在y=0處)，如圖7～圖9所示。由圖7～圖9可知，本文提出的改進模型可以獲得更接近于CONTACT的非線性剪切應力分布結果，從而進一步驗證了本文所提出的改進模型在提高接觸計算精度方面的有效性。

圖7 縱向剪切應力對比圖(y=0)

圖8 橫向剪切應力對比圖(y=0)

圖9 總剪切應力對比圖(y=0)

4.3 接觸力預測

本節(jié)旨在評估簡化方法預測接觸力的有效性。接觸力作為車輛動力學仿真的重要輸入量，其準確地進行預測對于車輛動力學仿真計算，尤其是車輛過曲線或道岔等特殊工況時有重要影響。

三種工況下滲透量與接觸力的關系曲線如圖10所示。在前兩種工況下，兩種簡化的接觸模型結果與CONTACT能夠較好地吻合。在最后一種工況下，盡管修正STRIPES模型顯著改善了STRIPES的計算精度，但仍存在進一步的改善空間。

圖10 不同橫移量下滲透量與接觸力關系的對比圖

純縱向蠕滑率變化時產生的縱向蠕滑力曲線，如圖11所示。修正STRIPES模型在前兩種工況可以獲得更接近CONTACT的計算結果，相比之下，原始STRIPES中蠕滑力的增加速度更緩慢。由于最后一種工況 Δy=6 mm中兩種接觸模型所得接觸斑形狀與CONTACT相比存在一定差異，兩種模型所得蠕滑力曲線的偏差較另外兩種工況偏大。

圖11 純縱向蠕滑工況不同橫移量時蠕滑力的對比圖

除純蠕滑情況外，還對純自旋條件下的蠕滑曲線進行了研究，如圖12所示。與純蠕滑情況類似，修正STRIPES模型結果更接近CONTACT計算的參考值。

圖12 純自旋工況不同橫移量下蠕滑力的對比圖

5 結論

簡化接觸模型STRIPES已廣泛應用于鐵路車輛軌道動力學和輪軌接觸計算，本文針對該模型提出了一些改進方法，并在修正STRIPES的模型中實現(xiàn)，該模型可與原模型保持相同的計算速度，并遠遠快于Kalker的CONTACT程序。通過與STRIPES模型和CONTACT程序計算結果進行詳細比較，可得出以下結論：

(1)使用變化的滲透量比例系數(shù)替代恒定比例系數(shù)這一改進對于改善法向接觸計算是必要且有效的，特別是對于整個接觸區(qū)域中存在橫向曲率變化的情況。

(2)本文提出了一種可靠并可廣泛用于其他非赫茲接觸模型的滲透迭代方法，提高了法向計算精度。

(3)在切向計算方面，改進STRIPES模型采用與FaStrip相結合的新型非橢圓方法，可顯著提高剪切應力分布和蠕滑力的計算精度，并獲得更加準確的蠕滑曲線。

(4)以上結果表明該方法具有更高的精確度，并可進一步應用于車輛軌道動力學和輪軌磨損計算，提高分析精度。

在下一步的研究中，本文將所提出的模型與車輛軌道動力學模型相結合，并繼續(xù)深入研究其對于磨耗輪軌型面的計算精度。