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

秦佳良，劉林芽，2

(1.華東交通大學 軌道交通基礎設施性能監測與保障國家重點實驗室，江西 南昌 330013；2.萍鄉學院 工程與管理學院，江西 萍鄉 337055)

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

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

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

1 車輛-軌道非線性耦合系統動力學方程

本文建立車輛-軌道非線性耦合系統分析模型時,由于車輛子系統和軌道子系統沿線路方向左右對稱,因此本文只取半結構展開研究。且本文僅考慮輪軌豎向動力效應,輪軌間考慮為Hertz非線性接觸。

1.1 車輛子系統動力學方程

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

車輛子系統的位移向量XV為

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

(1)

由Lagrange方程我們可以直接得到車輛子系統的動力學方程為

(2)

QV=-g{Mc0Mt0Mt0MwMwMwMw}T

(3)

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

(4)

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

(5)

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

1.2 軌道子系統動力學方程

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

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

(6)

(7)

最后利用有限元法的“對號入座”法則,根據板式無砟軌道單元的質量矩陣、阻尼矩陣和剛度矩陣,組集整個軌道子系統的質量矩陣、阻尼矩陣、剛度矩陣,可得到軌道子系統的動力學方程為

(8)

2 顯隱式求解算法計算流程

在用數值方法求解車輛-軌道非線性耦合系統分析模型的動力學方程時,因車輛子系統的質量矩陣為對角陣,數值計算時可以采用新型顯式積分法來求解車輛子系統的動力學方程;軌道子系統采用有限元法建立,其質量矩陣為非對角陣,數值計算時可以采用Newmark-β積分法來求解軌道子系統的動力學方程。分別建立車輛-軌道非線性耦合系統振動分析的分離迭代和分離同步兩種求解方法。分離迭代法是將車輛-軌道非線性耦合系統劃分為車輛和軌道兩個非時變的子系統模型,通過輪軌位移協調或輪軌力平衡關系將兩個子系統模型聯系起來,兩個子系統模型之間通過迭代計算來分別求解車輛子系統模型和軌道子系統模型的動力學方程。分離同步法其實也是將車輛-軌道非線性耦合系統分析模型劃分為車輛和軌道兩個非時變子系統模型,但車輛子系統模型和軌道子系統模型的動力學方程是同步求解的,不需要進行迭代求解。

2.1 顯隱式積分求解格式

(9)

(10)

式中:Δt為時間積分步長;ψ、φ分別為控制積分方法特性的兩個獨立參數,計算起步時只需令ψ=φ=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-β積分常數α和δ分別取0.25和0.5。

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

在開始求解時先假設車輛子系統和軌道子系統從靜止起步,即初始時刻t=0時,車輛子系統和軌道子系統的位移響應、速度響應和加速度響應均為零。在對時間步長進行循環計算時,假設在t+Δt時刻,兩個子系統動力學方程求解已進行k次迭代,現在考察第k+1次迭代。

當k為偶數時,令

(14)

并計算差值Δ1(i)為

(15)

當k為奇數時,令

(17)

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

Step6根據軌道子系統位移差值Δ進行收斂性判別,即

(18)

收斂條件為

(19)

Step7如果計算滿足收斂性,則可以轉入下一時間步長后繼續循環計算,直至整個時域T。如果計算不滿足收斂性,則令k=k+1,進入下一迭代步后繼續循環計算。

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

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

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

3 顯隱式求解算法驗證

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

圖3 本文計算結果

圖4 參考文獻[15]中計算結果

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

4 顯隱式求解算法收斂性對比分析

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

4.1 時間步長的確定

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

圖5 輪軌振動加速度的數值模擬結果

圖6 輪軌力的數值模擬結果

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

4.2 計算效率的對比分析

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

表1 不同時間步長耦合系統求解時間

由表1可知,采用顯隱式分離迭代法求解時,當時間步長取得較大時,采用Aitken加速法可以增強系統迭代求解的計算穩定性,但當時間步長取得越小,Aitken加速法的作用會越來越小。這是因為當時間步長較大時,耦合系統在下一時刻的平衡狀態會和上一時刻耦合系統的平衡狀態產生嚴重偏離的情況,這將導致耦合系統數值解的發散而使得計算結果無法收斂,或者耦合系統需要多次迭代才能使計算結果收斂,此時采用Aitken加速法可以減少迭代次數,可以起到增強迭代計算穩定性的作用。而當時間步長取得越小時,耦合系統在下一時刻的平衡狀態會和上一時刻耦合系統的平衡狀態越接近,這會使得耦合系統非常容易地再次達到其平衡位置,因此耦合系統的收斂速度會較快且Aitken加速法的作用會越來越小。

當采用相同時間步長時,顯隱式分離同步法的計算效率要比顯隱式分離迭代法更高。這是因為采用顯隱式分離迭代法求解時,每一時間步內車輛子系統和軌道子系統的動力學方程至少分別求解2次才能收斂,而采用顯隱式分離同步法求解時每一時間步內只需求解一次耦合系統的動力學方程,因此顯隱式分離同步法的計算效率會更高。但是顯隱式分離迭代法可以通過選取較大的時間步長來提高耦合系統的計算效率,在實際工程應用時可根據實際情況選用合適的求解方法。

5 結論

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

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

3)當時間步長較大時,采用Aitken加速法可以起到增強顯隱式分離迭代法的計算穩定性的作用。但隨著時間步長的減小,Aitken加速法的作用會越來越小。

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