車輛-軌道非線性耦合系統(tǒng)顯隱式求解算法的收斂性分析

秦佳良，劉林芽，2

(1.華東交通大學(xué) 軌道交通基礎(chǔ)設(shè)施性能監(jiān)測與保障國家重點(diǎn)實(shí)驗(yàn)室，江西 南昌 330013；2.萍鄉(xiāng)學(xué)院 工程與管理學(xué)院，江西 萍鄉(xiāng) 337055)

鐵路交通是國家發(fā)展的重大公共基礎(chǔ)設(shè)施,是現(xiàn)代化綜合交通運(yùn)輸體系的基礎(chǔ)骨干,對(duì)促進(jìn)經(jīng)濟(jì)社會(huì)的快速發(fā)展起到了十分顯著的作用。我國在進(jìn)入21世紀(jì)以來,在高速鐵路、重載鐵路、城市軌道交通等領(lǐng)域取得了世界矚目的成就。軌道交通通過輪軌相互作用來實(shí)現(xiàn)其功能,保持軌道良好的平順性,降低機(jī)車車輛對(duì)軌道和軌下基礎(chǔ)的動(dòng)力作用,減少養(yǎng)護(hù)維修的成本,以及新型機(jī)車車輛和新型軌道結(jié)構(gòu)的設(shè)計(jì)、制造等問題的解決都需要對(duì)車輛與軌道系統(tǒng)之間相互作用的動(dòng)力特性展開深入研究[1]。因此準(zhǔn)確、高效地求解車輛-軌道耦合系統(tǒng)的動(dòng)態(tài)響應(yīng)對(duì)保證軌道交通的可持續(xù)發(fā)展具有重大的理論和現(xiàn)實(shí)意義。

車輛-軌道耦合系統(tǒng)的動(dòng)力響應(yīng)模型求解一直是車輛動(dòng)力學(xué)和軌道動(dòng)力學(xué)研究中最關(guān)鍵的基礎(chǔ)問題[2]。近幾十年以來,國內(nèi)外眾多學(xué)者曾先后對(duì)車輛-軌道耦合作用引起的機(jī)車車輛、軌道結(jié)構(gòu)、路基、隧道、橋梁及大地振動(dòng)等方向展開了深入的研究[3-7],提出一系列繁簡各異的輪軌動(dòng)力學(xué)理論、算法、模型、程序及軟件,部分實(shí)現(xiàn)了工程應(yīng)用,有效地解決了在軌道交通可持續(xù)發(fā)展中面臨的一些瓶頸問題,也取得了一些較為顯著的研究成果[8-12]。在一般應(yīng)用情況下,車輛-軌道耦合系統(tǒng)的動(dòng)力學(xué)分析模型較復(fù)雜且動(dòng)力自由度數(shù)目也較多,采用解析法直接分析和求解其動(dòng)力學(xué)微分方程組會(huì)十分復(fù)雜且困難,因此目前一般都普遍采用直接數(shù)值積分法來求解。直接數(shù)值積分法分為顯式數(shù)值積分法和隱式數(shù)值積分法兩種不同的形式。隱式數(shù)值積分法一般具有穩(wěn)定性較好、計(jì)算的精度較高等特點(diǎn),普遍使用的隱式數(shù)值積分法一般包括Newmark-β法、Wilson-θ法、Park法和Hobolt法等。但隱式數(shù)值積分法每步積分均需求解大型代數(shù)方程組,這對(duì)計(jì)算規(guī)模大的工程問題帶來了巨大的挑戰(zhàn)。與隱式數(shù)值積分法不同,顯式數(shù)值積分法則一般具有計(jì)算過程較簡單、計(jì)算效率較高等一些優(yōu)點(diǎn),普遍使用的顯式數(shù)值積分法一般包括新型快速顯式積分法、四階Runge-Kutta法和中心差分法等。顯式數(shù)值積分法明顯的缺點(diǎn)就是其穩(wěn)定性或計(jì)算精度一般比隱式數(shù)值積分法略差?？梢?顯式數(shù)值積分與隱式數(shù)值積分的優(yōu)缺點(diǎn)是互補(bǔ)的。鑒于此,我們?cè)谇蠼廛囕v-軌道非線性耦合系統(tǒng)振動(dòng)方程時(shí),嘗試將顯式數(shù)值積分與隱式數(shù)值積分聯(lián)合使用,以期獲得比較好的綜合性能。

本文利用有限元法分別建立了車輛子系統(tǒng)和軌道子系統(tǒng)的動(dòng)力學(xué)分析模型,然后采用顯式數(shù)值積分法和隱式數(shù)值積分法分別求解車輛子系統(tǒng)和軌道子系統(tǒng)的動(dòng)力學(xué)方程,提出了車輛-軌道非線性耦合系統(tǒng)振動(dòng)分析的分離迭代和分離同步兩種求解方法,并對(duì)兩種方法進(jìn)行了算例驗(yàn)證。最后通過進(jìn)一步研究時(shí)間步長對(duì)輪軌耦合系統(tǒng)振動(dòng)響應(yīng)的影響,以及對(duì)比不同時(shí)間步長下車輛-軌道非線性耦合系統(tǒng)的數(shù)值計(jì)算效率,分析兩種算法在計(jì)算穩(wěn)定性、準(zhǔn)確性和數(shù)值編程求解效率等方面的差異,為求解車輛-軌道非線性耦合相互作用問題提供可供選擇的高效數(shù)值計(jì)算方法。

1 車輛-軌道非線性耦合系統(tǒng)動(dòng)力學(xué)方程

本文建立車輛-軌道非線性耦合系統(tǒng)分析模型時(shí),由于車輛子系統(tǒng)和軌道子系統(tǒng)沿線路方向左右對(duì)稱,因此本文只取半結(jié)構(gòu)展開研究。且本文僅考慮輪軌豎向動(dòng)力效應(yīng),輪軌間考慮為Hertz非線性接觸。

1.1 車輛子系統(tǒng)動(dòng)力學(xué)方程

車輛子系統(tǒng)可以考慮為一個(gè)10自由度的具有一、二系懸掛的整車模型,模型分別考慮車體、轉(zhuǎn)向架的沉浮和點(diǎn)頭運(yùn)動(dòng)以及輪對(duì)的沉浮運(yùn)動(dòng),具體情況見圖1。圖1中,Mc、Jc分別為車體的質(zhì)量、轉(zhuǎn)動(dòng)慣量;Mt、Jt分別為轉(zhuǎn)向架的質(zhì)量、轉(zhuǎn)動(dòng)慣量;Mwi(i=1,2)為第i個(gè)車輪的質(zhì)量;Ks1、Cs1分別為車輛的一系懸掛剛度、阻尼;Ks2、Cs2分別為車輛的二系懸掛剛度、阻尼;νc、νti(i=1,2)分別為車體、前后轉(zhuǎn)向架沉浮運(yùn)動(dòng)的豎向位移;θc、θti(i=1,2)分別為車體、前后轉(zhuǎn)向架點(diǎn)頭運(yùn)動(dòng)的角位移;νwi(i=1,2,3,4)為第i個(gè)車輪沉浮運(yùn)動(dòng)的豎向位移;l1、l2分別為車輛固定軸距之半、車輛定距之半。

車輛子系統(tǒng)的位移向量XV為

XV={vcθcvt1θt1vt2θt2vw1vw2vw3vw4}T

(1)

由Lagrange方程我們可以直接得到車輛子系統(tǒng)的動(dòng)力學(xué)方程為

(2)

QV=-g{Mc0Mt0Mt0MwMwMwMw}T

(3)

FVT={0 0 0 0 0 0F1F2F3F4}T

(4)

式中:g為重力加速度;Fi(i=1,2,3,4)為輪軌相互用力,可用Hertz非線性彈性接觸理論來計(jì)算,其表達(dá)式為

(5)

其中,νlci、ηi(i=1,2,3,4)分別為第i個(gè)輪軌接觸處的鋼軌位移、軌道的不平順幅值;G為輪軌接觸常數(shù)[5]。

1.2 軌道子系統(tǒng)動(dòng)力學(xué)方程

以路基上采用的CRTS Ⅱ型板式無砟軌道結(jié)構(gòu)為例,建立模擬板式無砟軌道子系統(tǒng)的動(dòng)力學(xué)分析模型。針對(duì)CRTS Ⅱ型板式無砟軌道子系統(tǒng)結(jié)構(gòu)形式的特點(diǎn),將軌道子系統(tǒng)沿線路縱向取半結(jié)構(gòu)進(jìn)行分析研究,并利用有限元法來模擬板式無砟軌道子系統(tǒng)。將沿線路縱向兩相鄰扣件之間的軌道結(jié)構(gòu)考慮為一個(gè)板式無砟軌道單元,單元長度記為l,單元模型見圖2。模型中,將鋼軌模擬為離散黏彈性點(diǎn)支承梁單元,Ky1、Cy1分別為扣件的彈性系數(shù)、阻尼系數(shù);將軌道板和底座板模擬為連續(xù)黏彈性支承梁單元,Ky2、Cy2分別為CA砂漿層的彈性系數(shù)、阻尼系數(shù),Ky3、Cy3分別為路基層的彈性系數(shù)、阻尼系數(shù)。圖2中,ν1、ν4分別為鋼軌兩端的豎向位移;θ1、θ4分別為鋼軌兩端的轉(zhuǎn)角;ν2、ν5分別為軌道板兩端的豎向位移;θ2、θ5分別為軌道板兩端的轉(zhuǎn)角;ν3、ν6分別為底座板兩端的豎向位移;θ3、θ6分別為底座板兩端的轉(zhuǎn)角。

圖2 板式無砟軌道單元模型

(6)

(7)

最后利用有限元法的“對(duì)號(hào)入座”法則,根據(jù)板式無砟軌道單元的質(zhì)量矩陣、阻尼矩陣和剛度矩陣,組集整個(gè)軌道子系統(tǒng)的質(zhì)量矩陣、阻尼矩陣、剛度矩陣,可得到軌道子系統(tǒng)的動(dòng)力學(xué)方程為

(8)

2 顯隱式求解算法計(jì)算流程

在用數(shù)值方法求解車輛-軌道非線性耦合系統(tǒng)分析模型的動(dòng)力學(xué)方程時(shí),因車輛子系統(tǒng)的質(zhì)量矩陣為對(duì)角陣,數(shù)值計(jì)算時(shí)可以采用新型顯式積分法來求解車輛子系統(tǒng)的動(dòng)力學(xué)方程;軌道子系統(tǒng)采用有限元法建立,其質(zhì)量矩陣為非對(duì)角陣,數(shù)值計(jì)算時(shí)可以采用Newmark-β積分法來求解軌道子系統(tǒng)的動(dòng)力學(xué)方程。分別建立車輛-軌道非線性耦合系統(tǒng)振動(dòng)分析的分離迭代和分離同步兩種求解方法。分離迭代法是將車輛-軌道非線性耦合系統(tǒng)劃分為車輛和軌道兩個(gè)非時(shí)變的子系統(tǒng)模型,通過輪軌位移協(xié)調(diào)或輪軌力平衡關(guān)系將兩個(gè)子系統(tǒng)模型聯(lián)系起來,兩個(gè)子系統(tǒng)模型之間通過迭代計(jì)算來分別求解車輛子系統(tǒng)模型和軌道子系統(tǒng)模型的動(dòng)力學(xué)方程。分離同步法其實(shí)也是將車輛-軌道非線性耦合系統(tǒng)分析模型劃分為車輛和軌道兩個(gè)非時(shí)變子系統(tǒng)模型,但車輛子系統(tǒng)模型和軌道子系統(tǒng)模型的動(dòng)力學(xué)方程是同步求解的,不需要進(jìn)行迭代求解。

2.1 顯隱式積分求解格式

(9)

(10)

式中:Δt為時(shí)間積分步長;ψ、φ分別為控制積分方法特性的兩個(gè)獨(dú)立參數(shù),計(jì)算起步時(shí)只需令ψ=φ=0。

(11)

(12)

(13)

式中:c0=1/(αΔt2);c1=δ/(αΔt);c2=1/(αΔt);c3=1/(2α)-1;c4=δ/α-1;c5=Δt(δ/α-2)/2;c6=Δt(1-δ);c7=δΔt;Newmark-β積分常數(shù)α和δ分別取0.25和0.5。

2.2 顯隱式分離迭代法求解流程

在開始求解時(shí)先假設(shè)車輛子系統(tǒng)和軌道子系統(tǒng)從靜止起步,即初始時(shí)刻t=0時(shí),車輛子系統(tǒng)和軌道子系統(tǒng)的位移響應(yīng)、速度響應(yīng)和加速度響應(yīng)均為零。在對(duì)時(shí)間步長進(jìn)行循環(huán)計(jì)算時(shí),假設(shè)在t+Δt時(shí)刻,兩個(gè)子系統(tǒng)動(dòng)力學(xué)方程求解已進(jìn)行k次迭代,現(xiàn)在考察第k+1次迭代。

當(dāng)k為偶數(shù)時(shí),令

(14)

并計(jì)算差值Δ1(i)為

(15)

當(dāng)k為奇數(shù)時(shí),令

(17)

注意到Aitken加速法是每迭代兩次才改善一次。

Step6根據(jù)軌道子系統(tǒng)位移差值Δ進(jìn)行收斂性判別,即

(18)

收斂條件為

(19)

Step7如果計(jì)算滿足收斂性,則可以轉(zhuǎn)入下一時(shí)間步長后繼續(xù)循環(huán)計(jì)算,直至整個(gè)時(shí)域T。如果計(jì)算不滿足收斂性,則令k=k+1,進(jìn)入下一迭代步后繼續(xù)循環(huán)計(jì)算。

為進(jìn)一步提高顯隱式分離迭代法的計(jì)算效率,參照文獻(xiàn)[14]中的改進(jìn)方法,根據(jù)式( 9 )、式(10)和式(11),可以在迭代步循環(huán)計(jì)算前,先計(jì)算車輛子系統(tǒng)的位移和速度,以及軌道子系統(tǒng)位移響應(yīng)在迭代步循環(huán)過程中的不變量,在迭代求解過程中直接調(diào)用即可,這樣可以減少計(jì)算工作量。

2.3 顯隱式分離同步法求解流程

與顯隱式分離迭代法一樣,顯隱式分離同步法在開始求解時(shí)車輛子系統(tǒng)和軌道子系統(tǒng)也從靜止起步,即假設(shè)初始時(shí)刻t=0時(shí),車輛子系統(tǒng)和軌道子系統(tǒng)的位移響應(yīng)、速度響應(yīng)和加速度響應(yīng)均為零?，F(xiàn)假設(shè)已經(jīng)求得t-Δt和t時(shí)刻車輛子系統(tǒng)和軌道子系統(tǒng)的位移響應(yīng)、速度響應(yīng)和加速度響應(yīng),現(xiàn)在需要求解t+Δt時(shí)刻兩個(gè)子系統(tǒng)的位移響應(yīng)、速度響應(yīng)和加速度響應(yīng)。主要計(jì)算步驟如下:

3 顯隱式求解算法驗(yàn)證

為驗(yàn)證本文提出的顯隱式分離迭代法和顯隱式分離同步法的正確性,與文獻(xiàn)[15]中計(jì)算的算例進(jìn)行比較。計(jì)算條件為假設(shè)列車編組為1動(dòng)+1拖的高速列車在博格板式無砟軌道上運(yùn)行,運(yùn)行速度為200 km/h,軌道高低不平順激振源考慮成波長為20 m、波幅為6 mm的周期性正弦函數(shù),其計(jì)算結(jié)果見圖3,文獻(xiàn)[15]的計(jì)算結(jié)果見圖4。

圖3 本文計(jì)算結(jié)果

圖4 參考文獻(xiàn)[15]中計(jì)算結(jié)果

由圖3和圖4的計(jì)算結(jié)果對(duì)比可知,鋼軌和博格板的位移響應(yīng)波形均符合物理概念,與文獻(xiàn)[15]計(jì)算得到的響應(yīng)幅值與變化規(guī)律基本一致,證明了本文提出的顯隱式分離迭代法和顯隱式分離同步法的正確性和可行性。

4 顯隱式求解算法收斂性對(duì)比分析

本文為對(duì)比分析顯隱式分離迭代法和顯隱式分離同步法的收斂性,以CRH3型高速動(dòng)車通過CRTSⅡ型板式無砟軌道為計(jì)算算例。其中車輛的運(yùn)行速度設(shè)置為300 km/h,板式無砟軌道單元的長度取0.65 m,在進(jìn)行數(shù)值計(jì)算時(shí)軌道子系統(tǒng)的長度設(shè)為360個(gè)軌道單元的長度。軌道高低不平順設(shè)置為車輛-軌道非線性耦合系統(tǒng)的激勵(lì)源,將其考慮成波長為12.5 m、波幅為3 mm的周期性正弦函數(shù)。車輛子系統(tǒng)和軌道子系統(tǒng)的結(jié)構(gòu)參數(shù)按照文獻(xiàn)[16]中的參數(shù)取值。

4.1 時(shí)間步長的確定

在應(yīng)用中最為重要的問題就是時(shí)間步長的確定,這直接影響到數(shù)值計(jì)算是否收斂以及計(jì)算效率。我們采用數(shù)值試驗(yàn)方法,對(duì)各種時(shí)間步長(最小時(shí)間步長為10-3ms)進(jìn)行了數(shù)值模擬,計(jì)算結(jié)果見圖5、圖6。

圖5 輪軌振動(dòng)加速度的數(shù)值模擬結(jié)果

圖6 輪軌力的數(shù)值模擬結(jié)果

由圖5和圖6可知,時(shí)間步長對(duì)輪對(duì)加速度和輪軌力的影響相對(duì)較小,而對(duì)鋼軌加速度的影響較大。顯隱式分離迭代法的臨界收斂時(shí)間步長為1 ms,超過1 ms后計(jì)算結(jié)果將不收斂,最大有效時(shí)間步長可取為0.5 ms。顯隱式分離同步法的臨界收斂時(shí)間步長和最大有效時(shí)間步長均為0.1 ms,超過0.1 ms后計(jì)算結(jié)果將不收斂。

4.2 計(jì)算效率的對(duì)比分析

為對(duì)比顯隱式分離迭代法和顯隱式分離同步法在計(jì)算效率上的差異,在計(jì)算機(jī)上運(yùn)行程序,計(jì)算時(shí)間步長分別取0.5、0.1、0.05、0.01 ms時(shí),考察兩種方法的計(jì)算耗時(shí)。不同時(shí)間步長下車輛-軌道非線性耦合系統(tǒng)的求解時(shí)間見表1。

表1 不同時(shí)間步長耦合系統(tǒng)求解時(shí)間

由表1可知,采用顯隱式分離迭代法求解時(shí),當(dāng)時(shí)間步長取得較大時(shí),采用Aitken加速法可以增強(qiáng)系統(tǒng)迭代求解的計(jì)算穩(wěn)定性,但當(dāng)時(shí)間步長取得越小,Aitken加速法的作用會(huì)越來越小。這是因?yàn)楫?dāng)時(shí)間步長較大時(shí),耦合系統(tǒng)在下一時(shí)刻的平衡狀態(tài)會(huì)和上一時(shí)刻耦合系統(tǒng)的平衡狀態(tài)產(chǎn)生嚴(yán)重偏離的情況,這將導(dǎo)致耦合系統(tǒng)數(shù)值解的發(fā)散而使得計(jì)算結(jié)果無法收斂,或者耦合系統(tǒng)需要多次迭代才能使計(jì)算結(jié)果收斂,此時(shí)采用Aitken加速法可以減少迭代次數(shù),可以起到增強(qiáng)迭代計(jì)算穩(wěn)定性的作用。而當(dāng)時(shí)間步長取得越小時(shí),耦合系統(tǒng)在下一時(shí)刻的平衡狀態(tài)會(huì)和上一時(shí)刻耦合系統(tǒng)的平衡狀態(tài)越接近,這會(huì)使得耦合系統(tǒng)非常容易地再次達(dá)到其平衡位置,因此耦合系統(tǒng)的收斂速度會(huì)較快且Aitken加速法的作用會(huì)越來越小。

當(dāng)采用相同時(shí)間步長時(shí),顯隱式分離同步法的計(jì)算效率要比顯隱式分離迭代法更高。這是因?yàn)椴捎蔑@隱式分離迭代法求解時(shí),每一時(shí)間步內(nèi)車輛子系統(tǒng)和軌道子系統(tǒng)的動(dòng)力學(xué)方程至少分別求解2次才能收斂,而采用顯隱式分離同步法求解時(shí)每一時(shí)間步內(nèi)只需求解一次耦合系統(tǒng)的動(dòng)力學(xué)方程,因此顯隱式分離同步法的計(jì)算效率會(huì)更高。但是顯隱式分離迭代法可以通過選取較大的時(shí)間步長來提高耦合系統(tǒng)的計(jì)算效率,在實(shí)際工程應(yīng)用時(shí)可根據(jù)實(shí)際情況選用合適的求解方法。

5 結(jié)論

1)提出了基于顯隱式積分格式的車輛-軌道非線性耦合系統(tǒng)分離迭代和分離同步兩種數(shù)值計(jì)算方法,并通過算例驗(yàn)證了兩種算法的正確性。

2)對(duì)車輛-軌道非線性耦合系統(tǒng)進(jìn)行數(shù)值分析時(shí),顯隱式分離迭代法和顯隱式分離同步法的最大有效時(shí)間步長分別是0.5、0.1 ms。

3)當(dāng)時(shí)間步長較大時(shí),采用Aitken加速法可以起到增強(qiáng)顯隱式分離迭代法的計(jì)算穩(wěn)定性的作用。但隨著時(shí)間步長的減小,Aitken加速法的作用會(huì)越來越小。

4)相同時(shí)間步長下顯隱式分離同步法的計(jì)算效率要比顯隱式分離迭代法高,但是可以采用較大的時(shí)間步長來提高顯隱式分離迭代法的計(jì)算效率。