輪軌蠕滑計(jì)算模型的比較研究

輪軌蠕滑計(jì)算模型的比較研究

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