許玉德，嚴(yán)道斌，孫小輝

（1. 同濟(jì)大學(xué)道路與交通工程教育部重點(diǎn)實(shí)驗(yàn)室，上海201804；2. 同濟(jì)大學(xué)上海市軌道交通結(jié)構(gòu)耐久與系統(tǒng)安全重點(diǎn)實(shí)驗(yàn)室，上海201804；3. 比亞迪汽車工業(yè)有限公司，廣東深圳518118）

列車的支承、轉(zhuǎn)向、加減速等行為都是由輪軌接觸區(qū)域的輪軌力完成的，輪軌接觸產(chǎn)生的應(yīng)力是導(dǎo)致輪軌滾動接觸疲勞和磨耗的重要因素，特別是在道岔區(qū)和曲線地段［1-2］。磨耗的發(fā)展導(dǎo)致輪軌接觸幾何關(guān)系惡化，進(jìn)而影響列車運(yùn)行的安全性和平穩(wěn)性。因此，準(zhǔn)確計算輪軌接觸應(yīng)力，對于了解輪軌作用機(jī)理，研究車輛-軌道耦合作用，評估列車運(yùn)行的安全性和平穩(wěn)性有著重要意義。

國內(nèi)外學(xué)者對輪軌接觸理論及計算模型開展了深入的研究。Hertz 理論［3］被最早運(yùn)用于輪軌接觸應(yīng)力計算，但該理論僅適用于平面或近似平面上的橢圓形、單點(diǎn)接觸。Kalker［4］提出了三維彈性體滾動接觸理論，基于此開發(fā)的Contact程序能夠較好地計算輪軌非橢圓形、多點(diǎn)接觸，是迄今為止使用最為廣泛的理論，但該理論使用了平面彈性半空間假設(shè)，因此不能較好地解決非平面接觸問題。Ayasse等［5］提出了基于虛位移原理的“條帶法”計算輪軌法向力以及利用Fastsim 程序［6-7］計算切向力，考慮接觸區(qū)域曲率對自旋的影響，但僅適用于固定曲率的接觸問題。Li［8］基于Kalker的三維彈性體滾動接觸理論，提出了利用四分之一空間研究輪軌接觸幾何的方法，開發(fā)的Wear程序適用于共形接觸的求解，但影響系數(shù)的計算存在精度不足的問題。Vollebregt 等［9］利用有限元方法計算了輪軌共形接觸的影響系數(shù)，解決了彈性半空間和四分之一空間假設(shè)的不足，并將該計算方法加入到Contact 程序中，但這使得Contact程序的計算效率有所下降。Zhu等［10］考慮了輪軌廓形局部曲率變化，并提出了一種改進(jìn)的半赫茲輪軌接觸計算模型。目前，有限元方法是計算輪軌接觸應(yīng)力最為準(zhǔn)確的方法，但計算效率過低制約了其用于批量計算的可能性［11］。

在道岔區(qū)或曲線地段，鋼軌磨耗通常較為嚴(yán)重，輪軌接觸面曲率變化較大，導(dǎo)致輪軌非平面接觸，進(jìn)而加劇磨耗，最終對鋼軌的使用壽命產(chǎn)生嚴(yán)重影響。準(zhǔn)確計算輪軌非平面接觸應(yīng)力，對預(yù)測磨耗、優(yōu)化廓形、延長輪軌使用壽命、保證列車平穩(wěn)通過［12-14］具有重要意義?；贙alker 的三維彈性體滾動接觸理論，結(jié)合輪軌非平面接觸幾何關(guān)系，對最小余能方程中的影響系數(shù)進(jìn)行了修正，提出了一種最小余能方程求解算法，實(shí)現(xiàn)了輪軌非平面接觸的法向應(yīng)力、切向應(yīng)力和黏著區(qū)-滑動區(qū)的計算。

1 理論方法

1.1 三維彈性體滾動接觸理論

如圖1 所示，在Kalker 的三維彈性體滾動接觸理論中，將可能接觸區(qū)域劃分為矩形單元網(wǎng)格，利用彈性半空間Boussinesq-Cerruti 公式計算得到該區(qū)域內(nèi)任一單元與其他所有單元的位移差，以法向余能最小為目標(biāo)，通過法向接觸算法Norm［4］求解實(shí)際接觸區(qū)域內(nèi)的法向應(yīng)力分布；以切向余能最小為目標(biāo)，通過切向接觸算法Tang［4］求解實(shí)際接觸區(qū)域內(nèi)的切向應(yīng)力分布以及黏著區(qū)-滑動區(qū)分布。

三維彈性體滾動接觸理論中最小余能方程為

圖1 可能接觸區(qū)域網(wǎng)格劃分Fig.1 Meshing of the possible contact region

式中：C1為余能；Ac為可能接觸區(qū)域；h為法向間隙；μn為法向位移；pn為法向接觸力；wτ為切向剛性滑移量；μτ為切向位移；μ′τ為前一時刻切向位移；pτ為切向接觸力；g 為滑動摩擦系數(shù)；P 為作用在接觸面上的豎向集中力。

對于平面接觸問題的求解，該理論認(rèn)為任一單元的位移與接觸面上作用力的對應(yīng)關(guān)系為

式中：A 為可能接觸區(qū)域各單元的影響系數(shù)；Ann為在法向方向作用單位力時產(chǎn)生的法向位移，即法向影響系數(shù)；Aττ為在切向方向作用單位力時產(chǎn)生的切向位移，即切向影響系數(shù)。

1.2 影響系數(shù)的修正

在三維彈性體滾動接觸理論中，法向影響系數(shù)Ann的計算采用的是彈性半空間的Boussinesq 解析解，切向影響系數(shù)Aττ的計算則采用的是Cerruti解析解。在輪軌非平面接觸問題中，彈性半空間假設(shè)不再成立，這就造成了Kalker的Contact程序無法解決輪軌非平面接觸問題，因此有必要對非平面接觸時的影響系數(shù)進(jìn)行修正。

參考文獻(xiàn)［15］的做法，如圖2所示，引入作用力F與曲面的法向角α用于影響系數(shù)修正，圖中符號含義、具體修正方法及步驟與文獻(xiàn)［16］完全相同，并在文獻(xiàn)［16］的基礎(chǔ)上完善了切向力作用時的影響系數(shù)修正公式。

圖2 影響系數(shù)修正Fig.2 Modification of the influence coefficient

根據(jù)文獻(xiàn)［15-16］的思路，當(dāng)可能接觸區(qū)域內(nèi)單元J 受到單位作用力時，可能接觸區(qū)域內(nèi)任一單元I產(chǎn)生的影響系數(shù)可表示為

式中：AIbJa為由Boussinesq-Cerruti 公式計算得到的影響系數(shù)，表示在單元J 受a 方向的單位作用力時，單元I 在b 方向的影響系數(shù)；為修正后的影響系數(shù)；a，b=1，2，3。為方便說明，將法向n、橫切向s、縱切向t分別用數(shù)字1、2、3表示，如即為，后文中的方向均按此方法表示。

1.3 影響系數(shù)矩陣的構(gòu)造及運(yùn)用

影響系數(shù)修正后，將可能接觸區(qū)域離散為由面積Δx1×Δx2的矩形構(gòu)成的網(wǎng)格，設(shè)橫切向網(wǎng)格數(shù)量為M，縱切向網(wǎng)格數(shù)量為N，共計MN 個矩形網(wǎng)格。將離散的二維矩形網(wǎng)格從1至MN進(jìn)行編號，在網(wǎng)格J（J=1，2，···，MN）的a（a=1，2，3）方向作用單位力時，記可能接觸區(qū)域內(nèi)所有網(wǎng)格的b（b=1，2，3）方向的影響系數(shù)

此時，可能接觸區(qū)域內(nèi)所有網(wǎng)格的影響系數(shù)可構(gòu)造為MN×MN維的影響系數(shù)矩陣，如下所示：

記沿著a（a=1，2，3）方向作用在可能接觸區(qū)域所有網(wǎng)格內(nèi)的接觸力

則所有網(wǎng)格在b（b=1，2，3）方向產(chǎn)生的位移可表示為

根據(jù)以上論述可知，通過影響系數(shù)矩陣構(gòu)造，對于任意輪軌接觸力的作用，利用式（3）～（7）就能計算離散后接觸面上所有單元的位移，進(jìn)一步可以計算該區(qū)域內(nèi)任一單元與其他所有單元的位移差。結(jié)合目標(biāo)函數(shù)，可以實(shí)現(xiàn)輪軌接觸區(qū)域法向應(yīng)力、切向應(yīng)力和黏著區(qū)-滑動區(qū)的計算。

2 算法實(shí)現(xiàn)

2.1 最小余能方程離散化

在計算輪軌非平面接觸時，不再以法向余能最小和切向余能最小為目標(biāo)函數(shù)，而是將總余能最小作為目標(biāo)函數(shù)，求解Kalker的最小余能方程。

將修正后的影響系數(shù)代入式（1）中，寫成離散形式并展開，所得到的即為輪軌非平面接觸條件下的最小余能方程，如下所示：

式中：hI為網(wǎng)格I（I=1，2，…，MN）處法向間隙為網(wǎng)格I 處a（a=1，2，3）方向位移；為網(wǎng)格I 處a（a=1，2，3）方向接觸力；為網(wǎng)格I處b（b=2，3）方向剛性滑移量為網(wǎng)格I 處b（b=2，3）方向前一時刻位移；αI為網(wǎng)格I處法向角。

此時，將所有變量表征為向量或矩陣的形式，則在輪軌非平面接觸條件下，最小余能方程的離散形式如 下所示：

式 中 ： h 為 法 向 間 隙 矩 陣 ， 即h=[h1h2… hMN]；w 為剛性滑移量矩陣，即w=[w1w2… wMN]。

2.2 最小余能方程求解分析

在式（9）中，p(1)、p(2)、p(3)3組接觸力需要求解，從目標(biāo)函數(shù)和約束條件來看，這是一個非線性約束規(guī)劃問題。首先，將式（9）中約束條件定義如下：

式中：?為自定義符號，表示相同維度的2 個矩陣之間相同位置元素對應(yīng)相乘，得到的是相同維度的新矩陣；cos α和sin α分別為法向角余弦值和正弦值矩陣 ， 即 cos α=[cos α1cos α2… cos αMN]，sin α=[sin α1sin α2… sin αMN]。

對于式（10），按照最優(yōu)解的Kuhn-Tucker 必要條件［17］，必然存在互補(bǔ)變量乘子≥0，≥0，∈R，使得下式同時成立：

式中：?為偏導(dǎo)符號；λ(a)為互補(bǔ)變量乘子矩陣，即

將式（9）中的C3分別對p(a)I求偏導(dǎo)，即?C3=，此時式（11）中的L1可全部展開，計算式如下所示：

式（11）、（12）中的6 個方程L11、L12、L13、L2、L3、L4即構(gòu)成了最優(yōu)解的Kuhn-Tucker 必要條件展開式，為非線性方程組。6 個方程中變量的數(shù)量包括：p(1)I、p(2)I、p(3)I、λ(1)I、λ(2)I各有MN 個，λ(3)I有1 個，總數(shù)量為（5MN+1）個。獨(dú)立方程的數(shù)量包括：L11、L12、L13各有MN個，L2、L3各有MN個，L4有1個，總數(shù)量為（5MN+1）個。變量與獨(dú)立方程的數(shù)量相同，因此可以進(jìn)行最優(yōu)化求解。

2.3 最小余能方程求解算法

直接采用牛頓-拉夫遜方法（N-R 法）［18］對6個方程進(jìn)行最優(yōu)求解是困難的，因?yàn)榉匠蘈2是一個三次方程，而N-R 法只是二次收斂的。對于類似的=0 互補(bǔ)松弛方程，和中必然有一個因子為0，因此可以通過判別，使其中一個因子為0，從而將方程降為二次方程，同時將方程簡化。

將整個計算區(qū)域分為接觸區(qū)域和非接觸區(qū)域。對于計算區(qū)域內(nèi)的任一網(wǎng)格I（I=1，2，…，MN），在接觸區(qū)域內(nèi)的所有法向接觸力必然有＞0，同時又要滿足互補(bǔ)松弛方程L3=0，因此在接觸區(qū)域內(nèi)對應(yīng)的=0。在非接觸區(qū)域內(nèi)，又有=0，則＞0。因此，當(dāng)假定了接觸區(qū)域后，就可以將互補(bǔ)松弛方程L3釋放掉。

對于接觸區(qū)域，又分為黏著區(qū)和滑動區(qū)。對于在接觸區(qū)域內(nèi)的任一網(wǎng)格I，在黏著區(qū)內(nèi)必然有，同時又要滿足互補(bǔ)松 弛 方 程，因此在黏著區(qū)內(nèi)對應(yīng)的=0。在滑動區(qū)內(nèi)，又有，則。通過 以上判定將L2降為二次方程，就能較好地使用N-R法。算法的基本步驟如下所示：

步驟1對可能接觸區(qū)域進(jìn)行網(wǎng)格劃分，輸入已知參數(shù)，以此計算影響系數(shù)、法向間隙等基本參數(shù)。

步驟2假定可能接觸區(qū)域全部為接觸區(qū)域，并設(shè)置為黏著區(qū)，此時的初始解均為0。

步驟3計算式（11）中的線性方程L1和L4，得到的初始解。

步驟4檢查接觸區(qū)域＞0是否成立。若是，則執(zhí)行步驟（5）；若不是，則將≤0的單元設(shè)定為非接觸區(qū)域，執(zhí)行步驟（8）。

步驟5檢查非接觸區(qū)域＞0 是否成立。若是，則執(zhí)行步驟6；若不是，則將≤0 的單元設(shè)定為接觸區(qū)域，執(zhí)行步驟8。

步 驟 6檢 查 黏 著 區(qū)是否成立。若是，則執(zhí)行步驟7；若不是，則將的單元設(shè)定為滑動區(qū)，同時設(shè)定切向接觸力

步驟7檢查滑動區(qū)＞0是否成立。若是，則當(dāng)前解即為計算結(jié)果，計算結(jié)束；若不是，則將≤0的單元設(shè)定為黏著區(qū)，執(zhí)行步驟8。

步驟8以上步驟中的作為初值T0，首先計算Tk（k=0，1，2，…）對應(yīng)的函數(shù)值L(Tk)和雅克比矩陣JL(Tk)，其次按照迭代公式ΔTk=-JL(Tk)-1L(Tk)計算迭代步長。若|ΔTk|≤ξ成立，其中ξ為精度條件，則執(zhí)行步驟3；若不成立，則更新變量使Tk+1=Tk+ΔTk，重新執(zhí)行步驟8，直到|ΔTk|≤ξ成立為止。

3 仿真驗(yàn)證

為驗(yàn)證所提算法的準(zhǔn)確性，以LM 型車輪和CHN75鋼軌為對象，利用所提算法計算輪軌軌距角非平面接觸時的法向應(yīng)力、切向應(yīng)力和黏著區(qū)-滑動區(qū)。同時，在有限元軟件中建立完全一致的工況并計算對應(yīng)的結(jié)果，通過對2 種算法計算結(jié)果的比較來判斷所提算法的準(zhǔn)確性。設(shè)置的工況及對應(yīng)的參數(shù)如表1 所示。需要說明的是，不考慮鋼軌截面縱向的變化。

表1 驗(yàn)證參數(shù)Tab.1 Parameters used for validation

如圖3 所示，在橫移量為+5.5 mm 條件下，車輪與鋼軌發(fā)生軌距角非平面接觸，該接觸為兩點(diǎn)接觸，計算結(jié)果可用于驗(yàn)證。

圖3 輪軌接觸二維斷面Fig.3 Two-dimensional cross section of wheel-rail contact

3.1 法向應(yīng)力

如圖4所示，從法向應(yīng)力計算結(jié)果來看，所提算法與有限元方法得到的結(jié)果是接近的。所提算法相比于有限元方法，在“高峰”處最大值差異為-2.6%，在“低峰”處最大值差異為-6.0%。

3.2 切向應(yīng)力

圖4 法向應(yīng)力結(jié)果Fig.4 Results of normal stress

如圖5所示，從切向應(yīng)力計算結(jié)果來看，所提算法與有限元方法得到的結(jié)果也是接近的。所提算法相比于有限元方法，在“高峰”處最大值差異為-2.0%，在“低峰”處最大值差異為-6.0%。

圖5 切向應(yīng)力結(jié)果Fig.5 Results of shear stress

3.3 黏著區(qū)-滑動區(qū)

如圖6 所示，從黏著區(qū)-滑動區(qū)分布的情況來看，所提算法相比于有限元方法，在“高峰”處的分布是比較接近的。值得注意的是，在“低峰”處有限元方法計算結(jié)果中出現(xiàn)了黏著區(qū)，但所提算法計算結(jié)果中并未出現(xiàn)。進(jìn)一步對黏著區(qū)-滑動區(qū)的面積進(jìn)行了對比，結(jié)果如表2 所示。所提算法相比于有限元方法，差異基本都控制在了7%以內(nèi)，這種差異還是可以接受的。

圖6 黏著區(qū)-滑動區(qū)分布Fig.6 Distribution of stick-slip region

表2 黏著區(qū)-滑動區(qū)面積Tab.2 Area of stick-slip region 單位：mm2

綜合上述對比，可以認(rèn)為所提算法能夠有效計算輪軌非平面接觸時的法向應(yīng)力、切向應(yīng)力和黏著區(qū)-滑動區(qū)。需要指出的是，由于有限元方法計算耗時過長，在計算本驗(yàn)證工況時超過1 h，而所提算法計算時長極短（不超過1 s），在計算效率上，所提算法更有優(yōu)勢，對車輛-軌道耦合動力學(xué)這類需要批量計算輪軌非平面接觸問題的求解也更為有利。

4 算例

利用所提算法研究鋼軌在不同磨耗狀態(tài)下的輪軌接觸特性。以我國某重載鐵路開行的27 t軸重列車為例進(jìn)行計算，選取通過總重約為60、90、200 MGT（1 MGT=106t）的磨耗鋼軌，取自同一曲線外軌測點(diǎn)。磨耗鋼軌型面如圖7所示，計算參數(shù)如表3所示。

圖7 磨耗鋼軌型面Fig.7 Worn rail profiles

表3 算例參數(shù)Tab.3 Parameters used for calculation

需要說明的是，第2、3次測量間隔時間較長，因此本算例僅是對不同磨耗狀態(tài)下的輪軌接觸特性變化趨勢做出判斷和分析。

4.1 法向應(yīng)力

如圖8 所示，在+13.0 mm 的橫移量條件下，隨著磨耗的增加，在60 MGT 時為典型的兩點(diǎn)接觸。雖然在90 MGT 時仍然屬于兩點(diǎn)接觸，但是由明顯的“雙峰”曲面向“單峰”曲面變化。到200 MGT時，接觸狀態(tài)過渡至共形接觸，應(yīng)力峰值出現(xiàn)在接觸中心。在此過程中，接觸斑形狀從60 MGT 時的縱向雙“橢圓形”，轉(zhuǎn)變至90 MGT時的橫向“葫蘆形”，直至200 MGT 時變成狹長的橫向“橢圓形”。在此過程中，最大法向應(yīng)力在逐漸減小，而接觸面積則逐漸增大。

圖8 法向應(yīng)力演變Fig.8 Evolution of normal stress

4.2 切向應(yīng)力

如圖9 所示，在+13.0 mm 橫移量的條件下，最大切向應(yīng)力在逐漸減小，值得關(guān)注的是，切向應(yīng)力峰值點(diǎn)的位置由60 MGT和90 MGT的雙接觸中心逐漸演變至200 MGT的多接觸中心。

圖9 切向應(yīng)力演變Fig.9 Evolution of shear stress

4.3 黏著區(qū)-滑動區(qū)

如圖10 所示，在+13.0 mm 橫移量條件下，黏著區(qū)-滑動區(qū)分布形狀出現(xiàn)了較為顯著的變化，黏著區(qū)的形狀由縱向“橢圓形”逐漸過渡至橫向“橢圓形”。表4列出了3種磨耗狀態(tài)下黏著區(qū)和滑動區(qū)的面積。可以看出，隨著總面積的增大，黏著區(qū)和滑動區(qū)的面積也逐漸增大。

圖10 黏著區(qū)-滑動區(qū)演變Fig.10 Evolution of stick-slip region

表4 算例的黏著區(qū)-滑動區(qū)面積Tab.4 Area of stick-slip region for calculation

5 結(jié)語

基于Kalker 的三維彈性體滾動接觸理論，結(jié)合輪軌非平面接觸幾何關(guān)系，提出了輪軌非平面接觸狀態(tài)下的影響系數(shù)計算公式，并構(gòu)造了影響系數(shù)矩陣，用于最小余能方程的離散化。將離散化最小余能方程的求解轉(zhuǎn)化為非線性規(guī)劃問題，根據(jù)對應(yīng)非線性方程組最優(yōu)解的Kuhn-Tucker 必要條件，提出了利用牛頓-拉夫遜方法求最小余能方程最優(yōu)解的算法，實(shí)現(xiàn)了輪軌非平面接觸狀態(tài)下法向應(yīng)力、切向應(yīng)力和黏著區(qū)-滑動區(qū)的計算。通過有限元方法驗(yàn)證了所提算法的準(zhǔn)確性，并利用該算法初步探索了鋼軌在不同磨耗狀態(tài)下的輪軌接觸特性。

該研究成果為精確、快速計算輪軌非平面接觸特性提供了新思路，未來將進(jìn)一步探索在輪軌磨耗預(yù)測、輪軌廓形優(yōu)化、輪軌使用壽命延長、列車運(yùn)行安全性和平穩(wěn)性評估等方面的適用性。