考慮基礎結構損傷的無砟軌道-車輛耦合動力模型及其求解

舒 瑤，蔣忠城，張 俊，張 波，劉國云，楊新文

(1. 西南交通大學牽引動力國家重點實驗室，四川，成都 610031；2. 大功率交流傳動電力機車系統集成國家重點實驗室，湖南，株洲 412001；3. 中車株洲電力機車有限公司，湖南，株洲 412001；4. 同濟大學道路與交通工程教育部重點實驗室，上海 201804)

無砟軌道結構損傷將會降低結構的承載能力，惡化輪軌關系，威脅運營安全，且軌道結構一旦出現傷損病害，不僅維修困難，而且維修費用昂貴。無砟軌道基礎結構主要為多孔性、多物相、極限應變小的水泥基準脆性材料，在車輛動荷載重復作用下，結構逐漸損傷劣化甚至斷裂失效。軌道結構損傷缺陷的動力演化及長時劣化，必然引起結構服役狀態與動態性能的持續改變[1 ? 3]，導致輪軌關系惡化，引起列車振動加劇，反過來又使得車輛系統對軌道的動力破壞作用增大，輪軌關系惡化與軌道結構損傷是一個長期的惡性循環過程。

國內外許多學者都曾圍繞混凝土結構損傷及其相關問題進行了研究[4 ? 12]，無砟軌道結構為典型的混凝土結構，為了研究無砟軌道結構的動力損傷問題，許多學者將混凝土結構分析理論引入無砟軌道傷損病害問題的研究，使無砟軌道結構的耐久性及傷損病害分析也大量采用數值計算同現代破壞力學理論相結合的分析方法。斷裂力學[3, 13 ? 14]、損傷力學[15 ? 16]、疲勞分析[16 ? 17]等現代破壞力學理論，以及界面力學[18 ? 19]理論和方法都被用于含損傷缺陷軌道的精細化分析。該類模型不再局限于線性的材料本構關系，而是引入斷裂力學、損傷力學等混凝土非線性力學行為的分析方法，建立精細化的數值模型，分析無砟軌道結構或其局部力學性能的變化。利用精細化數值模型進行動力仿真分析時，計算十分費時，且該類模型為了消除數值模型人工截斷邊界的影響，而將實際模型范圍延拓，往往導致計算量巨大。為深入、細致地研究軌道結構傷損的產生與動態演化，而不得不追本溯源地考慮輪軌之間的動態相互作用?，F有可以考慮軌道結構損傷缺陷發展演化的精細化數值模型因其計算量過大，難以實現車-軌耦合動力學分析，這也是目前很少見到同時考慮材料非線性及接觸非線性的車輛-軌道耦合動力分析的原因之一。

筆者曾在文獻[20 ? 21]中提出了一種能考慮軌道結構損傷發展的非線性軌道動力學模型，并利用擬譜法(Pseudo-spectral method，PSM)實現了模型的高效快速求解。為了探究無砟軌道結構損傷與車軌動態相互作用的相互影響機制，本文進一步基于該考慮軌道結構損傷的動力學模型，將車輛簡化為多剛體系統，且假定隨機性的軌面不平順以單節車長為周期重現，利用Hertz 非線性接觸實現車輛系統與含損傷無砟軌道系統的垂向傳力耦合，從而建立含損傷無砟軌道-車輛垂向耦合周期性動力模型；為了加速該同時包含材料非線性和接觸非線性的動力模型求解的收斂速度，采用隱式動力預測-校正算法和軌道-車輛系統交叉迭代的求解策略，實現了含損傷無砟軌道-車輛垂向耦合動力模型的隱式快速求解。

1 模型及運動方程

圖1 所示為含損傷無砟軌道-車輛垂向耦合系統，包含上部車輛子系統和下部含損傷無砟軌道子系統。將車輛簡化為包括車體、轉向架、輪對在內的多剛體系統，具有車體沉浮和點頭運動、前后轉向架沉浮和點頭運動、4 個輪對垂向運動等10 個自由度[22 ? 23]，單節高速車輛運動方程為：

圖1 含損傷無砟軌道-車輛垂向耦合系統Fig.1 Vertical coupling system of ballastless track and vehicle with damage

層狀的無砟軌道體系在豎向上具有極大的材料不均勻性和差異性，自上而下的鋼軌、軌道板、墊層和支承層的變形性能和損傷破壞模式差異巨大，在線彈性范圍內評估軌道結構的動力響應，忽視了不同結構材料的損傷特性和及其承載能力的差異，因而難以精確表征軌道系統的動力行為特性。無砟軌道的軌道板、支承層均為受拉易開裂的水泥基準脆性材料，研究表明，該類材料基本上都是帶缺陷或損傷服役，而且在列車荷載、溫度力等荷載或環境因素長期作用下，無砟軌道結構的損傷會不斷發展而影響其服役性能。真實材料的變形往往既有剛度退化，又有塑性變形存在，無砟軌道混凝土基礎結構在列車動荷載循環加卸載作用下動應變不大，一般忽略塑性變形，僅考慮彈性應變，軌道結構在服役過程中由于損傷累積而導致剛度降低的現象，則利用混凝土彈性損傷本構來描述。故從宏觀唯象的角度考慮混凝土受拉損傷的影響，采用《混凝土結構設計規范》(GB 50010?2010，2015 年版)[24]推薦的混凝土單軸受拉的彈性損傷本構模型，該彈性損傷本構的數學關系式可以表達為：

該彈性損傷本構應力-損傷-應變曲線如圖2 所示，無砟軌道軌道板和支承層的初始損傷場往往難以給定，本文采用的是方法是在荷載作用下，利用式(2)計算生成，軌道板C55 損傷因子d0為0.175 891，支承層C15 損傷因子d0為0.062 208。

圖2 彈性損傷本構應力-損傷-應變曲線Fig.2 Stress-damage-strain curves for elastic damage constitutive

如圖1 所示，截取周期長度為車長LC的一段無砟軌道結構，鋼軌、軌道板及支承層均采用歐拉梁模擬。鋼軌、軌道板和支承層的振動方程[20 ? 21]分別見式(3)～式(5)：

車輛運動方程式(1)與軌道系統振動方程通過輪軌接觸關系互相耦合，輪軌間Hertz 非線性接觸輪軌力pj(t)為[22 ? 23, 25]：

式中：相對位移δZ=Zwj(x,t)?[zr(xpj,t)+zr0(xpj,t)]，輪軌接觸常數G=3.86R?0.115×10?8m·N?2/3(磨耗型踏面車輪)；Zwj(x,t)為t 時刻第j 位輪對的垂向位移；zr(xpj,t)為t 時刻第j 位輪對下鋼軌的垂向位移；zr0(xpj,t)為t 時刻第j 位輪對下軌道的不平順。

如圖3 所示，假定隨機性的軌面不平順以單節車長為周期重現，進行車輛運動仿真計算時，先假設Car A 完全占據單個周期長度的軌道節段，如圖3(a)所示，隨著列車向前行進，鄰車Car B也逐漸進入該軌道節段，如圖3(b)～圖3(c)所示，直至鄰車Car B 完全占據該軌道節段，圖3(a)～圖3(d)恰好為單節車輛通過周期軌道節段的全過程。仿真過程中涉及到2 節車輛系統的平衡，利用模型的周期性進行簡化，當Car A 超出周期性軌道模型的邊界時，將鄰車Car B 前轉向架下的車輪所受到的輪軌力作為Car A 前轉向架下的車輪所受到的輪軌力，軌道系統的車輛荷載也需進行對應處理。

圖3 含損傷無砟軌道-車輛系統耦合仿真計算示意圖Fig.3 Schematic diagram of coupling simulation of ballastless track vehicle system with damage

2 模型求解

2.1 含損傷無砟軌道振動方程空間域離散

文獻[20 ? 21]應用周期性譜方法縮減模型求解規模，利用基函數的周期性處理模型周期性的邊界條件，并給出了含損傷無砟軌道振動方程的離散過程。譜方法(Spectral method)起源于Ritz-Galerkin 方法，以逼近性質良好、且便于利用快速變換來計算的正交多項式(Fourier 多項式、Chebyshev多項式、Legendre 多項式等)作為基函數，又分為Galerkin 譜方法、Tau 方法或擬譜方法(Pseudospectral method，又稱譜配點法)。譜方法利用正交基函數系將待求函數展開作為其某種意義上的逼近函數，然后通過人為截斷，把無限維問題簡化為有限維問題，將截斷的逼近函數展式代入原偏微分方程，利用Galerkin 方法將偏微分方程轉化為常微分方程，而后求解常微分方程的得到逼近函數的系數，從而求得原方程解。

鋼軌受到移動輪載作用，應用周期性Fourier-Galerkin 譜方法(簡記為FGSM)比較方便，將鋼軌位移寫成Fourier 級數，利用Galerkin 法將鋼軌振動偏微分方程離散成常微分方程組。

進而應用擬譜方法對軌道板和支承層的振動方程進行空間離散，令軌道板和支承層的振動方程在每一個計算配點滿足加權殘值為0，從而將軌道板和支承層的振動方程離散為常微分方程組。最后將離散后的常微分方程組整合為：

圖4 計算配點和應力點示意圖Fig.4 Schematic diagram of calculation collocation points and stress points

2.2 含損傷無砟軌道-車輛耦合動力模型的時間域求解

車輛-軌道耦合系統中，車輛系統被簡化為由質量-彈簧-阻尼集成的線性剛柔體系，其運動方程的求解采用常規多剛體運動學求解方法即可；而軌道系統被簡化為考慮材料損傷非線性的彈性地基上的疊合梁，其振動方程的求解往往需要采用增量迭代的方法，結構應力-應變的更新采用映射回退算法。而輪軌間Hertz 非線性接觸進一步引入了車輛-軌道耦合系統運動方程邊界條件的非線性，在每一時間步又必須采用迭代的方法求解。為了加速該同時包含材料非線性和接觸非線性的動力模型求解的收斂速度，采用隱式動力預測-校正算法和軌道-車輛系統交叉迭代的求解策略，以期實現含損傷無砟軌道-車輛垂向耦合動力模型的隱式快速求解。

2.2.1 隱式動力預測-校正算法

為了加速考慮損傷效應的無砟軌道振動方程的求解，在隱式求解方法[21]的基礎上，進一步引入預測-校正求解策略。動力方程的增量形式[21]可寫為：

由預估值式(10)計算得到不平衡的殘余力后，繼續利用殘余力計算得到廣義位移增量，并利用其對預估值進行校正，如此重復操作，直至達到前后兩個迭代子步廣義位移的計算結果之差滿足預設容差，此時認為迭代收斂而轉入下一時間步繼續計算。t+?t時步內迭代計算過程概要如下：

① 令迭代計數變量i=1，開始預估階段，令：

③ 用下式計算等效廣義剛度矩陣 K?：

⑤ 進入修正階段：

⑥ 如果第i 迭代步不是初始迭代步，判斷第i 步與第i?1 步兩個迭代子步軌道系統動位移的計算結果之差是否滿足預設容差，不滿足則令i=i+1，并轉到第②步，滿足則轉向⑦；

⑦ 利用式(16)得到的修正值作為t+?t時刻軌道的最終計算結果，并進一步更新軌道的應變及損傷等其他物理量，而后開始下一時間步的計算。

2.2.2 軌道-車輛系統隱式動力交叉迭代求解

為便于根據方程特性針對性地選用求解方法，分別對車輛系統及軌道系統交叉求解，如圖4所示，即在每一時間步，首先通過假設車輛系統、軌道系統的動位移，分別提取車輪位移、輪軌接觸點的軌面高低不平順及鋼軌位移，試算Hertz 非線性接觸輪軌力，進而將其輸入軌道系統求解軌道結構的動力響應而得到新的鋼軌位移，更新Hertz非線性接觸輪軌力后將其代入車輛系統求解車輛的動力響應而得到新的車輪位移，如此便完成了一個非線性接觸的循環迭代過程，判斷是否已經滿足迭代收斂條件，不滿足則再次更新輪軌力繼續迭代計算，反之則進入下一時間步的計算。

嚴格來說，考慮材料損傷非線性的軌道系統運動方程的求解，包含兩重意義的非線性問題的迭代，如圖5 所示，分別是軌道系統的平衡迭代和本構關系的迭代，其動力響應的計算需采用增量迭代的方法，另外車輛-軌道耦合的Hertz 非線性接觸也需要迭代，再加上時域求解固有的最外層時間步循環，故整個車輛-軌道耦合系統運動方程的時域數值求解的完整過程包含4 層循環，且其中3 層是非線性迭代。為簡化求解過程，嘗試將這三個層面的非線性迭代計算合并，不再需要在每一次試算輪軌力作用下的軌道系統都達到完全的迭代平衡，而是在Hertz 非線性接觸的迭代計算中不斷根據殘余力修正軌道系統的位移，以期達到材料本構、軌道系統及車輛-軌道耦合系統三個層面非線性迭代的平衡。如此處理雖然簡化了求解流程，但要達到迭代收斂往往會增加迭代步數。

圖5 每一時步迭代計算示意圖Fig.5 Schematic diagram of iterative calculation for each time step

為了加速迭代收斂，借鑒Newton-Raphson 法的思想，將隱式動力預測-校正算法嵌入輪軌非線性接觸的循環迭代過程中，提出一種同時求解材料非線性和接觸非線性的預測-校正迭代算法。在時刻t+?t，車輛系統和軌道系統動力方程的典型形式分別為：

同理，t+?t 時刻車輛系統的位移預估值為：

進而利用軌道、車輛系統的位移預估值，試算Hertz 非線性接觸輪軌力，然后展開t+?t時步各迭代子步的計算，直到達到前后兩個迭代子步軌道系統動位移的計算結果之差滿足預設容差，此時認為迭代平衡而轉入下一時間步繼續計算，t+?t時步內迭代計算過程概要如下：

1) 令迭代計數變量i=1，開始預估階段，令軌道系統的位移、速度和加速度為式(11)，令車輛系統的位移為式(19)；

3) 按照第2.2.1 節所述預測-校正算法中t+?t 時步迭代計算過程②～⑤，計算軌道系統的位移，更新過程變量及Hertz 非線性接觸輪軌力，進而將其輸入車輛系統求解車輛的動力響應而得到新的車輪位移；

4) 如果第i 步不是初始迭代步，判斷第i 步與第i?1 步兩個迭代子步軌道系統動位移的計算結果之差是否滿足預設容差，不滿足則令i=i+1并轉到第2)步，滿足則轉向第5)步；

5) 利用步驟4)得到的輪軌力生成軌道系統廣義荷載向量，再按②～⑤修正軌道系統的廣義位移、速度和加速度，以便得到與車輛系統位移同步更新的軌道系統位移，將其作為t+?t時刻軌道系統的最終計算結果，并進一步更新軌道系統的應變及損傷等其他物理量，而后開始下一時間步的計算。

軌道-車輛系統隱式動力交叉迭代求解流程見圖6。

3 模型驗證與討論

3.1 預測-校正算法的驗證與討論

為了驗證預測-校正算法的正確性，將其計算結果與隱式動力N-R 法[21]的計算結果進行比較，時間步長Δt =1×10?4s。圖7 給出了預校正算法與隱式動力N-R 法軌道板和支承層位移計算結果的對比，圖7 中x 表示縱向坐標，t 表示時間，ΔZs表示軌道板垂向位移絕對差，ΔZb表示支承層垂向位移絕對差。由圖可以看出，兩種方法的軌道板和支承層位移計算結果的絕對差均非常小，可以認為兩種方法的計算結果一致，預校正算法的位移計算結果準確可靠。

圖6 交叉迭代求解示意圖Fig.6 Schematic diagram of cross iteration solution

進一步提取預測-校正算法與隱式動力N-R 法軌道板和支承層最終損傷計算結果進行對比，如圖8 所示，圖中x 表示縱向坐標，t 表示時間，Δds表示計算的軌道板損傷因子絕對差，Δdb表示計算的支承層損傷因子絕對差。兩種方法的軌道板和支承層最終損傷計算結果的絕對差也都非常小，可以認為兩種方法最終損傷的計算結果一致，進一步驗證了預校正算法計算結果的正確性，也相互應證了兩種計算方法的可靠性。

3.2 含損傷無砟軌道-車輛垂向耦合動力模型的驗證與討論

以和諧號高速動車組CRH3 為例[25]，進行含損傷無砟軌道-車輛耦合動力分析，CRH3 高速動車組一般為8 節編組，最高運營速度達到350 km/h，仿真所用單節車輛系統參數見表1，軌道系統參數見表2。

圖7 軌道板和支承層位移計算結果對比Fig.7 Comparison of displacement calculation results of slab and bed

軌面高低不平順采用中國高速鐵路無砟軌道軌面高低不平順頻譜[27]，計算所采用的高低不平順樣本見圖9，假定隨機性的軌面不平順以單節車長為周期重現時，周期性模型偏向于考慮波長小于單節車長的不平順成分，如需考慮波長成分更為豐富的軌面不平順，可以將周期性模型的周期長度拓展為整數倍單節車輛長度。

為驗證含損傷無砟軌道-車輛垂向耦合動力模型計算的正確性，特將無砟軌道系統退化到彈性無損階段，采用FGSM 法[20 ? 21]計算了鋼軌為離散點支承的無砟軌道-車輛系統的動力響應，簡記為FGSM-e，也采用PSM 法[20 ? 21]計算了鋼軌為連續均勻支承的彈性無損階段無砟軌道-車輛系統的動力響應，簡記為PSM-e；而將含損傷無砟軌道-車輛垂向耦合動力模型簡記為PSM-d，三者均采用Hertz 非線性接觸模擬車輛-軌道耦合系統的動力相互作用。進而分別提取輪軌力F、車體加速度及軌道結構位移等物理量進行對比，從而驗證含損傷無砟軌道-車輛垂向耦合動力模型計算的正確性。圖10 分別給出了FGSM-e、PSM-e 和PSM-d 三種模型計算結果的對比，圖10 中x 表示縱向坐標，t 表示時間，F 表示輪軌力，由圖可以看出，FGSM-e、PSM-e 和PSM-d 三種模型輪軌力F、車體加速度及軌道結構位移等物理量的計算結果基本一致，相互應證了各自模型計算結果的正確性。

圖8 軌道板和支承層損傷計算結果對比Fig.8 Comparison of damage factors calculation results of slab and bed

表1 車輛系統參數Table 1 Parameters of vehicle system

表2 軌道系統參數Table 2 Parameters of track system

表3 給出了FGSM-e、PSM-e 和PSM-d 三種模型計算耗時對比，由表可知，對于彈性階段的無砟軌道-車輛垂向耦合動力模型的求解，PSM 法求解速度最快，僅耗時67.1 s，而FGSM 法耗時100.8 s；考慮損傷效應后，利用PSM 法求解含損傷無砟軌道-車輛垂向耦合動力模型需耗時125.8 s，每一時間步長內最大迭代11 次收斂，整個仿真時間內每一時間步長平均迭代6 次，計算仍然非常高效，方便應用于多參數、多工況的仿真計算。

圖10 計算結果對比Fig.10 Comparison of calculation results

表3 計算耗時對比Table 3 Calculation time comparison

圖11 軌道板縱剖面損傷分布云圖Fig.11 Damage distribution of longitudinal section of track slab

圖12 軌道板上、下表面應變與損傷因子變化曲線Fig.12 Time history curves of strain and damage factor on slab surface

圖11 給出了單節車輛通過后軌道板縱剖面損傷分布云圖，圖中x 表示縱向坐標，z 表示垂向坐標，以軌道板中心線為坐標0 位置，由該圖可以看出損傷更多主要集中在軌道板的上、下表面。進一步提取軌道板跨中上、下表面應變與損傷因子，圖12 給出了軌道板跨中上、下表面應變與損傷因子時程曲線，t 表示時間，由圖可以看出，在車輛動荷載作用下，相對于初始的損傷因子0.175 891，含初始損傷的軌道板的損傷因子會有所增加，在當前的車輛動荷載水平下，軌道板在列車動載作用下彎曲時，如圖12(b)所示，軌道板底部最大拉應變出現在車輛輪軸作用時刻，最大壓應變出現在相鄰兩節車輛轉向架作用的中間時刻，由最大負彎矩引起；如圖12(a)所示，軌道板頂部最大壓應變出現在車輛輪軸作用時刻，最大拉應變出現在相鄰兩節車輛轉向架作用的中間時刻，也由最大負彎矩引起，并導致了軌道板頂部受拉損傷的進一步增加，故在進行軌道結構設計時，負彎矩對軌道板頂部受拉損傷的影響不容忽視。

4 結論

本文運用車輛-軌道耦合動力學理論，利用混凝土彈性損傷本構考慮無砟軌道結構的損傷效應，將車輛簡化為多剛體系統，建立考慮軌道結構損傷效應的無砟軌道-車輛垂向耦合周期性動力模型，并采用隱式動力預測-校正算法和軌道-車輛系統交叉迭代的求解策略，實現了同時包含材料非線性及接觸非線性的含損傷無砟軌道-車輛垂向耦合動力模型的隱式快速求解。所建模型不僅計算結果正確，而且具有非常高的求解效率，方便應用于多參數、多工況的仿真計算，可用于分析和研究無砟軌道結構損傷與車軌動態相互作用的相互影響規律。