輪軌蠕滑計算模型的比較研究

輪軌蠕滑計算模型的比較研究

作為車輛動力學(xué)的基礎(chǔ)理論，f.c.chart首先提出了滑動的概念，并給出了平行同質(zhì)圓柱之間的縱向滑動力和縱向滑動率之間的非線性關(guān)系。1964年K.L.Johnson和P.J.Vermeulen對Carter的工作進行了拓展,提出了接觸面是橢圓、蠕滑率任意的縱、橫向蠕滑系數(shù)計算公式,但仍沒有考慮輪踏面的錐度。1967年J.J.Kalker運用三維彈性體理論,使輪軌蠕滑理論提升到了實用階段,提出了三維(縱向、橫向和回旋)蠕滑理論,蠕滑力與蠕滑率的關(guān)系是非線性的。其后編制了非線性的DUVOROL計算機程序和簡化的FASTSIM計算機程序,以及線性理論的計算公式,在80年代至90年代又推出了采用精確彈性理論的計算機程序CONTACT和USETAB數(shù)表。Z.Y.Shen和J.K.Hedrick等人又根據(jù)以上理論在1983年提出了實用的計算機公式化近似算法。近年來又提出了復(fù)橢圓的非赫茲接觸理論和輪緣接觸的磨合接觸。輪軌蠕滑力與蠕滑率之間的關(guān)系是很復(fù)雜的,它不僅與蠕滑率有關(guān),還與輪軌之間的接觸斑形狀、接觸壓力的分布、輪軌表面狀態(tài)及輪軌材質(zhì)等因素有關(guān)。因而計算蠕滑力時若要全面考慮這些因素,在工程應(yīng)用上幾乎是不可能的。又由于某些參量的離散性很大,對車輛動力學(xué)性能的影響又不敏感,所以精確計算輪軌蠕滑力從工程應(yīng)用角度考慮也是沒有必要的。所以以上各種算法都基于不同的假設(shè),適用于不同的場合。而且這些算法產(chǎn)生于不同的時代,受到了計算技術(shù)的影響。隨著PC機的迅猛發(fā)展,PC機無論在速度上,還是在內(nèi)存容量上都有了本質(zhì)的改變。如以前制約DOS程序的數(shù)組尋址限制問題在現(xiàn)在的Windows程序上已不再存在,計算速度與最初的PC286相比已提高了幾個數(shù)量級,內(nèi)存原來為640KB,現(xiàn)在可以選配128MB以上。因而有必要根據(jù)現(xiàn)有的計算條件和工程應(yīng)用的水平對蠕滑力模型和算法進行分析和評估,以更好地應(yīng)用于機車車輛的動力學(xué)分析計算。1輪軌歷史模型模擬兩個彈性滾動體在正壓力作用下,在相互接觸處的切平面上能產(chǎn)生切向力,在該切平面的法線方向能產(chǎn)生回轉(zhuǎn)力矩。這些力和力矩的大小主要與上下接觸面之間的蠕滑率等有關(guān),蠕滑率分為切向蠕滑率和回轉(zhuǎn)蠕滑率,切向蠕滑率在兩垂直的縱橫坐標(biāo)軸上的投影分別為縱向蠕滑率和橫向蠕滑率。切向力在同一坐標(biāo)軸上的對應(yīng)分量稱為縱向蠕滑力和橫向蠕滑力,如圖1所示。Fx、Fy和Mz分別為縱向蠕滑力、橫向蠕滑力和回旋蠕滑力矩。無量剛化后的各向蠕滑率定義為:vx=Vx′-VxVvy=Vy′-VyVωz=Ωz′-ΩzVvx=Vx′?VxVvy=Vy′?VyVωz=Ωz′?ΩzV蠕滑力除了與這些蠕滑率有關(guān)以外還與接觸區(qū)的形狀、正壓力的分布、接觸區(qū)的表面潤滑狀態(tài)、材質(zhì)特性參數(shù)等有關(guān)。它們之間的關(guān)系一般是非線性的。尋求合理的數(shù)學(xué)模型,并以工程允許的精度來反映它們之間的關(guān)系是輪軌蠕滑模型建模的關(guān)鍵,即在給定蠕滑率和其它參數(shù)及邊界條件下計算獲得具有一定精度的蠕滑力和蠕滑力矩?，F(xiàn)在鐵路輪軌接觸的蠕滑模型主要考慮單點接觸,并假定勻質(zhì)橢球體,即接觸區(qū)的形狀為橢圓。如圖1所示,橢圓長軸之半為a,短軸之半為b。由于解決問題的性質(zhì)不同和發(fā)展的漸進過程,常用的模型在形式上可大致分為公式化、復(fù)雜程序化和表格化三種。根據(jù)模型的細化程度又可以分為線性的、非線性的和準(zhǔn)非線性的。1.1旋轉(zhuǎn)活性公式1958年Johnson將Carter的兩維理論延伸到兩個滾動球體的三維工況,其中包含有縱向蠕滑和橫向蠕滑,但沒有自旋蠕滑。1964年Johnson與Vermeulen又將光滑的半空間理論引入研究沒有自旋蠕滑的純?nèi)浠r。其中提出了對合成蠕滑力的修正公式,但沒有計及回旋蠕滑的影響。1967年Kalker提出了計算蠕滑力的線性理論和計算公式,在小蠕滑率下依此計算有一定的精度,且計算方便。對于較大蠕滑率情況下的蠕滑力計算,則結(jié)合Kalker的線性理論公式和J-V的修正公式,得到了“縮減因子公式法”(又稱沈氏定理),它應(yīng)歸類于公式化的準(zhǔn)非線性法。首先用Kalker的線性公式計算線性蠕滑力Fx和Fy:Fx=-G(ab)C11vx,Fy=-G(ab)C22vy-G(ab)3/2C23ωz合成的切向蠕滑力F′r為:F′r=(Fx2+Fy2)12F′r=(Fx2+Fy2)12然后修正合成切向蠕滑力Fr:Fr={μΝ(F′rμΝ)-13(F′rμΝ)2+127(F′rμΝ)3(F′r≤3μΝ時)μΝ(F′r＞3μΝ時)Fr={μN(F′rμN)?13(F′rμN)2+127(F′rμN)3(F′r≤3μN時)μN(F′r＞3μN時)因Kalker的橫向蠕滑力計算中包含回轉(zhuǎn)蠕滑的影響,所以在修正公式中也包含了回轉(zhuǎn)蠕滑力矩的成分。但從公式的本身看是不精確的?；剞D(zhuǎn)蠕滑力矩Mz的計算公式為:Mz=G(ab)3/2C23vy-G(ab)2C33ωz式中:a、b分別為接觸橢圓的長短半軸長;G為剛度模量;μ為摩擦系數(shù);Ν為正壓力;C11,C22,C23,C33為蠕滑系數(shù);vx,vy,ωz分別為縱向、橫向和自旋蠕滑率。1.2基于高速算法的輪軌求解1978年Kalker根據(jù)三維精確理論編制了DUVOROL程序,但算法復(fù)雜,計算很慢。1982年編制了CONTACT程序,雖然增強了計算功能,但在計算速度上并沒有大的加快。同期Kalker提出基于新簡化理論的快速算法程序FASTSIM,為尋求速度,作了較大的簡化,但在模型中有錯誤。Kalker指出由于輪軌表面的狀態(tài)參數(shù)和摩擦系數(shù)有很大離散度,對輪軌間的切向力計算的精度要求不必太高。這是在當(dāng)時的計算機條件下提出的實用途徑,因而在80年代FASTSIM程序在曲線通過的計算中得到了廣泛應(yīng)用。1.3計算參數(shù)的壓縮自適應(yīng)識別插值表格法最初在1980年被英國Derby鐵路研究機構(gòu)所采用。隨著計算機內(nèi)存容量的大幅度提高和相應(yīng)軟件編譯器的發(fā)展,Kalker在1995年為了兼顧計算速度和精度又再次推出了插值表格法(USETAB),比英國Derby采用的表格方法精度更高。USETAB由4.5兆的數(shù)表TABCON文件和一些輔助程序組成,具體結(jié)構(gòu)如圖2所示。TABCON文件內(nèi)的數(shù)據(jù)由CONTACT程序計算產(chǎn)生,三列數(shù)據(jù)分別是作用在輪子上的規(guī)格化蠕滑力和蠕滑力矩,是規(guī)格化的蠕滑率和a/b的函數(shù)。柏松比為0.28。材質(zhì)的彈性模量G為82000N/mm2。規(guī)格化的蠕滑率為:vx′=-G(ab)c11vx3FΝvx′=?G(ab)c11vx3FNvy′=-G(ab)c22vy3FΝvy′=?G(ab)c22vy3FNωz′=-G(ab)32c23ωzFΝωz′=?G(ab)32c23ωzFN規(guī)格化的蠕滑力為:Fx′=FxFΝFx′=FxFNFy′=FyFΝFy′=FyFNΜz′=Μz√abFΝMz′=Mzab√FN該方法在1996年由Kalker引入我國,在所有快速算法中是精度最高的。目前在我國尚未推廣應(yīng)用。要將4.5兆的數(shù)據(jù)讀入內(nèi)存,采用基于DOS的640kB的內(nèi)存的Fortran編譯器是不可能實現(xiàn)的,在微機上可能的途徑是用NDPFortran,Fortran90或FortranPowerStation。另一方法是將該程序翻譯成其它支持大數(shù)組的語言,如VC++等。這里對USETAB作一簡單介紹。Kalker羅列了9個與計算蠕滑力有關(guān)的變量:vx、vy、ωz、a0/b0、c(√a0*b0)、G、si(柏松比)、μ和正壓力N。形成數(shù)表時將G、c、Nμ設(shè)置為常數(shù)1。將柏松比si設(shè)置為0.28,得到四個獨立的變量:a/b、vx、vy和ωz。將其中三個變量規(guī)格化為:vx′=-Gc2C11vx/(3Nμ)vy′=-Gc2C22vy/(3Nμ)ωz′=-Gc3C33ωz/(Nμ)整個表格根據(jù)a/b的取值平均分為1/7<a/b≤1和1≤a/b<7兩個部分,每個部分又根據(jù)s與a/b(在區(qū)間上1≤a/b<7上取b/a)的乘積取整后的值平均分成1～7七個單元;再將每個單元按照vx′的取值各分為0<vx′≤1和vx′>1兩個部分,然后每個部分再按照s與vx′(vx′>1時取vx′=1/vx′)的乘積取整后的值平均分成0～7八個單元;再將每個單元按照vy′的取值各分為0<vy′≤1和vy′<1兩個部分,然后每個部分再按照s與vy′(vy′>1時取vy′=1/vy′)的乘積取整后的值平均分成0～7八個單元;再將每個單元按照ωz′的取值各分為0<ωz′≤1、ωz′>1、0>ωz′>-1、和ωz′≤-1四個部分,然后每個部分再按照s與ωz′(ωz′≥1時取ωz′=1/ωz′;0>ωz′>-1時取ωz′=-ωz′;ωz′≤-1時取ωz′=-1/ωz′)的乘積取整后的值平均分成0～7八個單元。當(dāng)vx′<0時,有Fx=-Fx(|vx′|),當(dāng)vy′<0時,也有Fy=-Fy(|vy′|)。因而共有:2×7×2×8×2×8×4×8×3個數(shù)據(jù)。如表1所示。當(dāng)a/b大于7或小于1/7時,由于難以用線性插值法,仍用改進后的FASTSIM程序計算。2對于工況參量的影響為便于討論,以下分別用“J-V”、“FASTSIM”和“USETAB”符號來代表以上三種算法。首先是目標(biāo)參照對象問題,雖然人們在輪軌蠕滑方面做出了大量的工作,但不同的試驗,因條件不一,結(jié)果相差很大。理論模型也有多種,在學(xué)術(shù)屆Kalker1982年推出的CONTACT程序被公認為較細化的模型。而USETAB的數(shù)據(jù)表源于CONTACT程序,它們之間的誤差為1.5%。因而這里以USETAB作為相對參照對象。第二個問題是提取信息問題,由于比較的主要數(shù)據(jù)有縱向蠕滑力、橫向蠕滑力和回旋蠕滑力矩,而這些數(shù)據(jù)與多個參量有關(guān),故是一個多維空間問題。這里將這些參量分為三大類,第一類是主動參量,包括橫向蠕滑率和縱向蠕滑率;第二類是工況參量,包括接觸面積、形狀、摩擦系數(shù)、正壓力等;第三類是回旋蠕滑率。橫向蠕滑率和縱向蠕滑率是可以同時連續(xù)變化的參量,比較時必須同時考慮。工況類參量是固定的,比較時可以作為條件參數(shù)處理?；匦浠实拇笮Q定于接觸區(qū)的錐度和輪徑,實際上是可以連續(xù)變化的,但對于給定的車輪踏面,由于其外形的錐度變化范圍和接觸位置幾乎是固定的,可以將它處理成三個區(qū)間,如在踏面上取較小的值,在輪緣上取較大的值,輪緣根部取中等值。經(jīng)過這樣的處理,橫向蠕滑力和縱向蠕滑力的比較問題就成為在不同條件和工況下的兩維問題(橫向蠕滑率和縱向蠕滑率)。回旋蠕滑力矩的比較問題也可以處理成不同工況下的兩維問題(橫向蠕滑率和縱向蠕滑率)。關(guān)于回旋蠕滑力矩與縱向蠕滑率的關(guān)系,由于對稱性,可以看成是與合成切向蠕滑率的關(guān)系,即可以用單橫向蠕滑率來代替。至于取值,可以用與USETAB的絕對誤差和相對誤差的最大值。關(guān)于參量取值的范圍問題,這里根據(jù)實際可能出現(xiàn)的范圍選取。橫向蠕滑率和縱向蠕滑率,從零至完全飽和為止,一般為10%;回旋蠕滑率可以根據(jù)踏面錐度1/Δ和輪對實際滾動圓半徑R求得:ωz=sin(arctg(1/Δ)R1/Δ的變化范圍一般為:1/40～2.75(取輪緣角為70°)。對應(yīng)的Δ為0.36～40,這里取Δ為0.35～40。對于半徑R,考慮到不同大小的干線機車車輛車輪和磨耗后的大小變化,R的可能變化范圍為:360～625mm。若以360作為實際輸入?yún)?shù),將R的變化合并到Δ的變化當(dāng)中,則修正后的Δ的變化范圍為0.35～70。當(dāng)計算橫向、縱向蠕滑力時,取較小自旋蠕滑率,即Δ取為60。關(guān)于工況類參數(shù),這里根據(jù)通常出現(xiàn)的情況,選擇了以下6種組合,由3種蠕滑水平(全蠕滑、半蠕滑水平和1/4蠕滑水平)和2種接觸斑長短軸之比組成。2.1種活性水難以融合小旋流體的活性表2為小自旋蠕滑率情況(即Δ為60)下,3種蠕滑水平時的縱向和橫向蠕滑力的差異比較。計算條件是:輪對滾動圓半徑為525mm,車輪踏面橫斷面外形的半徑為無窮大,鋼軌頭部的半徑為300mm,接觸斑的長短軸之半分別為8.224mm和5.757mm。表3為小自旋蠕滑率情況(即Δ為60)下,三種蠕滑水平時的縱向和橫向蠕滑力的差異比較。計算條件是:輪對滾動圓半徑為525mm,車輪踏面橫斷面外形的半徑為450mm,鋼軌頭部的半徑為300mm,接觸斑的長短軸之半分別為8.224mm和5.757mm。2.2回轉(zhuǎn)園滑率值的計算縱向蠕滑率取為零,對應(yīng)于不同的橫向蠕滑率,當(dāng)回旋蠕滑率變化時,用USETAB、FASTSIM、J-V三種方法計算所得回旋蠕滑力矩數(shù)值的比較結(jié)果如表4所示。3輪軌探索的橫向動態(tài)曲線通過USETAB插值表格法、FASTSIM算法和飽和因子公式算法的比較表明,當(dāng)輪軌接觸點位于踏面上,即回旋蠕滑率較小時,不同算法之間的差異變化比較大。由FASTSIM算法求得