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

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

1 研究背景

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

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

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

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

2 STRIPES算法

2.1 法向接觸計算

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

(1)

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

x——滾動方向的坐標(biāo)；

y——橫向方向的坐標(biāo)。

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

(2)

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

V——材料的泊松比；

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

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

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

(3)

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

(4)

式中：A0——第一個接觸點(diǎn)處的縱向曲率；

B0——第一個接觸點(diǎn)處的橫向曲率；

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

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

(5)

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

接觸斑內(nèi)任意單元的應(yīng)力為：

(6)

2.2 切向接觸計算

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

圖1 STRIPES模型示意圖

(7)

(8)

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

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

υy——橫向蠕滑率；

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

(9)

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

G——剪切模量(MPa)。

假設(shè)牽引邊界遵循拋物線分布：

(10)

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

(11)

否則，該點(diǎn)位于滑動區(qū)，此時剪切應(yīng)力為：

(12)

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

3 STRIPES模型修正

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

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

(13)

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

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

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

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

(14)

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

δn=λδn-1

(15)

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

(16)

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

3.3 與FaStrip相結(jié)合的非橢圓方法

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

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

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

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

(17)

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

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

(18)

4 計算結(jié)果分析

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

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

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

4.1 法向接觸

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

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

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

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

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

4.2 切向接觸結(jié)果

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

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

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

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

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

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

4.3 接觸力預(yù)測

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

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

圖10 不同橫移量下滲透量與接觸力關(guān)系的對比圖

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

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

除純?nèi)浠闆r外，還對純自旋條件下的蠕滑曲線進(jìn)行了研究，如圖12所示。與純?nèi)浠闆r類似，修正STRIPES模型結(jié)果更接近CONTACT計算的參考值。

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

5 結(jié)論

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

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

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

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

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

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