基于DEM-MFBD方法的有砟軌道路基不均勻沉降影響分析

余翠英 ,雷紅博 ,羅文俊 ,馬斌

（1.華東交通大學(xué) a.理學(xué)院；b.土木建筑學(xué)院，南昌 330013；2.中鐵二局集團有限公司，成都 610013）

在鐵路運輸高速化、重載化背景下，有砟軌道因其經(jīng)濟性及適應(yīng)性較強，在世界范圍內(nèi)廣泛鋪設(shè)。由于地質(zhì)條件復(fù)雜及其基床動應(yīng)力較大[1]，容易發(fā)生路基不均勻沉降及軌枕空吊等病害[2]。路基不均勻沉降已成為影響有砟軌道結(jié)構(gòu)性能和耐久性的關(guān)鍵問題[3]。

傳統(tǒng)試驗方法成本高、可重復(fù)性差，難以觀測道床細觀力學(xué)行為[4]，當今大多數(shù)學(xué)者采用數(shù)值模擬方法開展軌道沉降受力變形及動力響應(yīng)研究[5]，尤其是離散元法，作為一種反映顆粒受力與運動狀態(tài)的數(shù)值方法，可準確表達散粒體道床的力學(xué)特性，被廣泛應(yīng)用于有砟軌道相關(guān)研究[6]。但是，單純離散元方法難以準確描述軌道結(jié)構(gòu)的連續(xù)體特征，部分學(xué)者采用離散元—有限元等耦合方法嘗試開展有砟軌道性能研究[7]。離散元—多柔性體動力學(xué)(DEM-MFBD)耦合方法作為一種可以有效反映離散體和連續(xù)體之間的接觸及其各自的物理力學(xué)特性的仿真方法，目前鮮有應(yīng)用于有砟軌道的相關(guān)研究。

針對有砟軌道上部結(jié)構(gòu)連續(xù)性與道床散粒體特征，筆者采用離散元—多柔性體動力學(xué)耦合方法，建立有砟軌道結(jié)構(gòu)沉降損傷模型，研究路基不均勻沉降對有砟軌道結(jié)構(gòu)整體受力變形的影響，以期為路基沉降控制和有砟軌道結(jié)構(gòu)綜合養(yǎng)護、維修提供理論依據(jù)和技術(shù)支持。

1 離散元與多柔性體動力學(xué)耦合方法

表1 為不同數(shù)值模擬方法的特性及適用范圍的對比。通過對比表1 中多種耦合方法的特性，采用離散元—多柔性體動力學(xué)耦合方法探討有砟軌道路基沉降大變形問題。

表1 不同數(shù)值模擬方法的對比Table 1 Comparison of numerical simulation methods

離散元接觸模型采用Hertz-Mindlin 非滑動模型[8]，該模型適于有砟軌道道砟非黏性顆粒的計算。由圖1 可見，道砟顆粒簇由若干個球單元鑲嵌而成，道砟顆粒簇之間的接觸即轉(zhuǎn)化為若干個球單元間的兩兩接觸，接觸單元間的接觸力由單元間的相互重疊量及相對速度決定。

圖1 Hertz-Mindlin 非滑動接觸模型Fig.1 Hertz-Mindlin (no slip) contact model

離散元—多柔性體動力學(xué)耦合求解中，墻單元網(wǎng)格面上的接觸力往往不在網(wǎng)格節(jié)點上，且每個三角形網(wǎng)格上可能有多個接觸點，需采用形函數(shù)插值法，結(jié)合虛功原理，將顆粒作用在耦合面上的接觸力轉(zhuǎn)換為等效節(jié)點力，其整體坐標與局部坐標轉(zhuǎn)換公式[9]為

式中：Ttrans3×3為轉(zhuǎn)換矩陣，則耦合邊界上等效接觸力計算公式[9]為

2 模型的建立

圖2 為DEM-MFBD 耦合過程流程圖，有砟軌道結(jié)構(gòu)建模過程分為離散元道床建模和多柔性體動力學(xué)上部結(jié)構(gòu)建模，通過軌枕墻單元作為耦合過渡邊界，進行軌道上部結(jié)構(gòu)和離散元道床耦合分析。

圖2 離散元?多柔性體動力學(xué)模型耦合流程圖Fig.2 Flowchart of coupling DEM-MFBD model

2.1 有砟軌道道床離散元模型

為準確模擬道床內(nèi)摩擦力和抗剪強度等力學(xué)行為，應(yīng)考慮道砟的真實幾何特征與粒徑級配[10-11]。選取華東交通大學(xué)軌道交通實驗基地有砟線路的特級道砟為樣本，利用三維掃描技術(shù)獲取的不規(guī)則道砟外形輪廓作為道砟顆粒外形邊界，填充外形邊界，建立道砟顆粒簇。如圖3 所示，根據(jù)《鐵路碎石道砟》(TB/T 2140—2008)[11]規(guī)范要求，首先剔除道砟針狀、片狀等異常樣本；其次，結(jié)合文獻[10,12]的經(jīng)驗并多次優(yōu)化計算，最終選取10 種道砟顆粒樣本，需要注意的是，樣本需借助NX 修改粒徑比例，以滿足顆粒級配要求，故選取樣本僅考慮其外形特征影響。其中，每個道砟顆粒平均由9.23 個大小不同的球單元構(gòu)成，用以平衡計算精度及計算效率。

圖3 道砟顆粒樣本建模過程Fig.3 The modeling process of ballast-particle sample

基于離散元的方法目前存在計算耗時高、占用內(nèi)存大等局限性，加上路基不均勻沉降問題對模型縱向長度的要求[13]，建立全尺度軌道模型難以同時滿足計算效率與精度的雙重要求，為此，需進行軌道橫斷面上的簡化。參考既有研究[6]并綜合考慮計算精度及計算效率，最終確定軌枕和道床橫斷面均縮尺為0.2 m 寬的2.5 維簡化模型，如圖4 所示。離散元道床模型長為6 m、寬為0.2 m、厚度為0.33 m，模型中共9 558 個道砟顆粒簇，考慮為特級道砟級配，如圖5 所示。

圖4 離散元道床模型Fig.4 The discrete element model of ballast bed

圖5 顆粒級配曲線Fig.5 Particle gradation curve

2.2 有砟軌道上部結(jié)構(gòu)多柔性體動力學(xué)模型

有砟軌道上部結(jié)構(gòu)采用多柔性體動力學(xué)建模（見圖6），將鋼軌、軌枕考慮為連續(xù)體，鋼軌視為柔性體，軌枕視為剛體。選用60 kg/m 鋼軌，長12.5 m，為實體單元；Ⅲ型彈條扣件，剛度為6×107N/m，采用bushing 單元模擬；Ⅲ型混凝土軌枕，軌枕間距為0.6 m，將軌枕作為耦合邊界導(dǎo)出為墻單元。為保證計算效率，僅將路基沉降區(qū)段6 m 范圍的道床進行離散元建模?？紤]到有砟道床對上部結(jié)構(gòu)的支承作用及不承受拉力的特點，且模型試驗中離散元區(qū)域外道床受力較小，考慮為無張力彈性，忽略其塑性變形。

圖6 離散元?多柔性體動力學(xué)耦合模型Fig.6 The coupling model by DEM-MFBD

模型簡化關(guān)鍵在于能否準確反映道床應(yīng)力狀態(tài)。假定列車荷載及上部結(jié)構(gòu)自重由軌枕傳遞至軌枕底面壓力為均勻分布，基于軌枕寬度不變，軌枕長度簡化為0.2 m，為保證道床頂面所受平均應(yīng)力不變，將軌枕底面壓力等比例縮減。相對于列車荷載，鋼軌和扣件自重影響較小，故僅考慮軌枕底面壓力受列車荷載單一因素影響，對輸入荷載進行等效折減。

式中：F為輸入總荷載；P為靜荷載；Fmax為動荷載幅值；k為縮減系數(shù)；f為輸入荷載頻率。

鋼軌外形尺寸與實際情況一致，由材料力學(xué)撓曲線近似微分方程可知，鋼軌抗彎剛度與彈性模量成正比，外荷載等效折減后應(yīng)對鋼軌進行剛度折減，以準確反映軌道的抗彎特性。

式中：E為鋼軌輸入彈性模量；E′為設(shè)計彈性模量。

3 參數(shù)標定及模型驗證

3.1 模型參數(shù)及標定

校準離散元模型參數(shù)的方法通?？煞譃閮煞N：1）基于微觀視角的直接測量法，直接在顆粒或接觸層面測量參數(shù)值；2）基于宏觀視角的方法，參數(shù)值在宏觀角度上反映顆粒堆積體的整體力學(xué)行為。顯然，直接測量法難以反映所有相關(guān)的物理效應(yīng)，在計算效率及精度的限制下，很難完全恢復(fù)顆粒的不規(guī)則形狀，且顆粒間的相互咬合也會被簡化。在宏觀方法中，不同參數(shù)組合可以得到相同的模擬結(jié)果，但參數(shù)的真實物理意義不明確[14]。為應(yīng)對參數(shù)組合的模糊性，需通過多種標準化試驗進行驗證。結(jié)合文獻[10,12]的參數(shù)選擇，在多次試算后確定模型參數(shù)，見表2。

表2 模型參數(shù)Table 2 Model parameters

為驗證計算參數(shù)反映道砟宏觀剪切特性的正確性，首先，參考文獻[10]的方法進行三軸壓縮試驗的數(shù)值模擬。如圖7 所示，三軸壓縮試驗樣本尺寸為0.4 m×0.2 m×0.2 m，以上述離散元道砟模型建立樣本，密度為1 410 kg/m3，孔隙率為45.8%。對模型施加50 kPa 的圍壓，將蓋板以緩慢的勻速運動下壓至應(yīng)變?yōu)?6.8%，記錄蓋板與道砟間的接觸力，以計算偏應(yīng)力及應(yīng)力比；進一步地，參照文獻[12]模擬多級剪切試驗，模型包括上下兩個箱體，在箱體內(nèi)部裝滿上述離散元道砟模型，樣本尺寸為0.3 m×0.3 m×0.2 m，密度為1 448 kg/m3，孔隙率為46.3%，對上方箱體向下階段性施加16.5、31.5、56.5、81.5、106.5 kPa 的壓強，并以0.01 m/s 的速度勻速平移下方箱體，記錄上方箱體側(cè)壁壓強，以計算切應(yīng)力。

圖7 數(shù)值樣本Fig.7 Numerical sample

水平應(yīng)力σ2=σ3，參考文獻[10]中的方法，通過平均有效應(yīng)力p'及差應(yīng)力q計算內(nèi)摩擦角?。

式中：σ1和σ2分別為軸向應(yīng)力和水平應(yīng)力。

圖8 所示為三軸壓縮試驗的差應(yīng)力及應(yīng)力比與豎向應(yīng)變的關(guān)系，圖9 所示為多級剪切試驗中的剪應(yīng)力變化，數(shù)值模擬計算抗壓試驗中內(nèi)摩擦角約為41°，多級剪切試驗中內(nèi)摩擦角為50.5°?？梢姡瑪?shù)值模擬結(jié)果與相關(guān)試驗[10,12]擬合較好，本文結(jié)果略大于文獻[10,12]結(jié)果，考慮到樣本級配的差異及本文道砟顆粒不可破碎，誤差在可接受范圍內(nèi)，說明選取的本構(gòu)模型可以還原道床的宏觀特性。

圖8 三軸壓縮試驗數(shù)值樣本的偏應(yīng)力及應(yīng)力比Fig.8 The differential stress and stress ratio of the numerical sample by triaxial compression test

圖9 多級剪切試驗數(shù)值樣本的應(yīng)力—應(yīng)變行為Fig.9 The stress-displacement behavior of numerical sample by multi-stage shear test

3.2 模型驗證

為驗證模型表達道床應(yīng)力狀態(tài)的正確性，基于提出的荷載折減方法，并考慮輪軌力為枕上壓力的3 倍[15]，對道砟箱試驗[16]工況進行模擬，將模型計算結(jié)果與試驗結(jié)果[6,16]進行比較。圖10 顯示了列車荷載沿軌道縱向5 根軌枕分散傳遞至道床（紅線表示接觸力走向，紅色越明顯表示接觸力越大），中間3根軌枕承擔(dān)75%以上的列車荷載，與既有研究結(jié)論一致[6]；模型測得最大接觸力為4 287.3 N，文獻[16]中最大接觸力為4 233 N，誤差約為1.3%，表明模型能夠反映真實道床應(yīng)力狀態(tài)。

圖10 道床接觸力鏈Fig.10 Contact force chain of ballast bed

為進一步驗證模型反映軌道結(jié)構(gòu)變形特性的正確性，參照同濟大學(xué)開展的有砟軌道路基不均勻沉降實尺模型試驗[4]，對鋼軌中心相距1.8 m 的兩點施加正弦荷載以模擬單一轉(zhuǎn)向架的動力作用。荷載分為靜載與動載兩部分，其中靜載為45 kN，動載幅值為49 kN。路基不均勻沉降為余弦型，路基沉降波長為2.4 m，沉降幅值為25 mm。通過試算并參考文獻[17]，最終確定荷載頻率為5 Hz，加載40 個完整周期，可得到鋼軌中心點隨加載過程的沉降曲線及基本穩(wěn)定狀態(tài)下的鋼軌沉降變形曲線。

由圖11 可以看出，鋼軌頂面中心點的豎直沉降量嚴格地隨著周期性循環(huán)荷載的作用而產(chǎn)生周期性變化，并產(chǎn)生彈性變形；由于鋼軌累計沉降量不斷增大，加載初期鋼軌沉降量較大，其增長趨勢隨循環(huán)加載次數(shù)的增加逐漸趨于平緩，與相關(guān)研究中對道床初期沉降與道床塑性變形的結(jié)論一致[18]。

圖11 軌面中心點隨循環(huán)加載過程的沉降曲線Fig.11 Settlement curve of central point of rail with loading process

如圖12 所示，將數(shù)值模擬得到鋼軌沉降變形與文獻[4]對比，僅兩數(shù)據(jù)點相對誤差較大，為11%和11.8%，其余數(shù)據(jù)點相對誤差范圍為0.2%～6.9%。數(shù)值模擬中鋼軌最大沉降量為7.47 mm，略小于文獻[4]中的結(jié)果，沉降曲線趨勢性擬合較好。考慮到試驗的離散性和趨勢性，模擬結(jié)果合理，數(shù)值模型可行，能夠反映道床真實應(yīng)力狀態(tài)、路基不均勻沉降下軌道結(jié)構(gòu)變形特性及沉降產(chǎn)生的歷程及機理。

圖12 鋼軌沉降變形曲線與文獻[4]的對比Fig.12 Comparison between the rail settlement deformation curve and literature [4]

4 計算結(jié)果分析

4.1 沉降波長對軌道受力變形的影響

為研究路基不均勻沉降波長下軌道的受力變形規(guī)律，鑒于文獻[4]足尺試驗工況中沉降波長均小于3 m，工況設(shè)定考慮波長范圍為5 根軌枕間距，故設(shè)定路基沉降幅值為10 mm，沉降波長分別為0.6、1.2、1.8、2.4 m 等4 種工況，加載50 個周期，荷載頻率為5 Hz，分析各路基沉降波長下軌道的結(jié)構(gòu)受力變形特征。

如圖13 所示，軌枕沉降與鋼軌沉降變形大致相同。路基沉降波長經(jīng)道床、軌枕傳遞到鋼軌表面會有一定程度的擴大。測得路基—鋼軌沉降波長擴散角θ為30.3°～34.4°，略大于文獻[4]試驗結(jié)果，結(jié)果表明，路基—軌面沉降波長擴散角受路基沉降波長、幅值變化影響較小，由于道床的塑性形變，路基沉降幅值及沉降面積傳遞至鋼軌表面均有減小。

圖13 沉降幅值為10 mm，波長為2.4 m 時軌道的整體變形Fig.13 The overall deformation of track structure with settlement amplitude10 mm and wavelength 2.4 m

如圖14 所示，路基沉降幅值為10 mm，波長由0.6 m 增大至2.4 m，鋼軌不平順幅值由1.82 mm 增大至6.34 mm，增大了約2.5 倍。同時，鋼軌兩端的上拱減少，鋼軌整體變形呈下沉趨勢。當路基沉降波長為0.6 m 時，鋼軌最大沉降量為?1.79 mm，而鋼軌上拱最高達到1.81 mm，由上拱導(dǎo)致的不平順較為嚴重。當路基產(chǎn)生小波長沉降時，鋼軌上拱導(dǎo)致的不平順將影響列車運營的舒適性。

圖14 不同路基沉降波長下鋼軌變形Fig.14 Rail deformation under different subgrade settlement wavelength

由圖15 可知，在路基沉降區(qū)相鄰軌枕處，道床應(yīng)力集中現(xiàn)象明顯，隨著路基沉降波長的增加，應(yīng)力集中位置由7 號軌枕逐漸外擴至8 號軌枕處，最大接觸力整體呈上升趨勢。正常情況下，最大接觸力為5 841.5 N，位于7 號軌枕下方；沉降波長為0.6、2.4 m 時，最大接觸力分別為7 227.8、10 405.0 N，分別上升了23.7%及78.1%。圖15 中數(shù)字為該位置軌枕與道床間接觸力，紅色突出顯示為其中最大的接觸力。

圖15 不同路基沉降波長下道床接觸力鏈Fig.15 Contact force chain of ballast bed under different subgrade settlement wavelengths

4.2 沉降幅值對軌道受力變形的影響

為研究路基不均勻沉降幅值對軌道受力變形的影響規(guī)律，基于工務(wù)維修規(guī)則和易操作性，同時借鑒文獻[4]的統(tǒng)計和現(xiàn)場調(diào)研，設(shè)定路基沉降波長為1.8 m，沉降幅值為10、15、20 mm 三種工況，加載情況同上，分析不同路基沉降幅值下軌道結(jié)構(gòu)受力變形特征。

如圖16 所示，隨著路基沉降幅值由10 mm 增大至15、20 mm，鋼軌不平順幅值分別增大了16%和29%。鋼軌沉降波長略有擴大，波長擴散角θ隨之減小，而端部上拱差別不大。

圖16 不同路基沉降幅值下的鋼軌變形Fig.16 Rail deformation under different subgrade settlement amplitudes

如圖17所示，隨著路基沉降幅值的增大，7號軌枕與道床的接觸力從7 679.3 N 逐漸降低至3 693.3 N，同時，8 號軌枕與道床接觸力從8 202.5 N 升高至10 827.5 N，相較于正常線路，上升了85.4%，應(yīng)力集中位置逐漸外擴（由7 號軌枕向8 號軌枕遷移）。

圖17 不同路基沉降幅值下的道床接觸力鏈Fig.17 Contact force chain of ballast bed under different subgrade settlement amplitude

4.3 路基沉降對軌枕空吊的影響規(guī)律

不同沉降波長及幅值的組合數(shù)值模擬結(jié)果如表3 所示。由表3 可知，軌枕空吊情況可由鋼軌沉降面積與路基沉降面積之比(S1/S0)劃分，當S1/S0>1.03 時，均沒有發(fā)現(xiàn)軌枕空吊；而當S1/S0≤1.03時，均產(chǎn)生不同程度的軌枕空吊。說明通過軌面沉降面積與路基沉降面積的比值反映軌枕空吊情況具有合理性和可行性，但對于該指標反映發(fā)生軌枕空吊的臨界值及影響機理仍需進一步研究。

表3 不同工況計算結(jié)果Table 3 The calculation results of different working condition

路基沉降波長為1.8 m 時，隨著沉降幅值由10 mm 增大至15、20 mm，S1/S0由1.22 下降至0.95和0.83，并產(chǎn)生了不同程度的軌枕空吊。沉降幅值為15 mm 時，軌枕在列車荷載作用下與道床接觸，而隨著列車荷載的減小，逐漸與道床分離，即為非完全空吊；沉降幅值繼續(xù)增大至20 mm 時，在列車荷載作用下，軌枕始終不與道床接觸，即為完全空吊。當路基沉降幅值為15 mm 時，沉降波長由0.6 mm 增大至2.4 mm，軌枕與道床間隙由0 mm 增大至4.3 mm，軌枕空吊加劇。

在路基沉降波長為0.6～2.4 m 的條件下，沉降幅值為10 mm 時，均未發(fā)生軌枕空吊；沉降幅值為15 mm 時，則發(fā)生了軌枕空吊。可見，路基沉降波長介于0.6～2.4 m 之間時，軌枕空吊臨界值的路基沉降幅值應(yīng)介于10～15 mm 之間，建議將路基沉降幅值控制在10 mm 以內(nèi)。

5 結(jié)論

1）基于DEM-MFBD 耦合方法，提出一種2.5維簡化模型，該模型可以正確反映有砟道床的宏觀力學(xué)行為，可應(yīng)用于有砟軌道沉降變形機理分析。

2）鋼軌沉降面積與路基沉降面積的比值可以反映軌枕空吊情況；當路基沉降波長范圍小于4 根軌枕間距時，軌枕空吊的沉降幅值臨界值介于10～15 mm 之間；建議沉降幅值控制在10 mm 以內(nèi)，避免軌枕空吊。

3）路基沉降波長和幅值的增加均導(dǎo)致軌道不平順明顯增大，并加劇軌枕空吊現(xiàn)象，建議加強路基沉降監(jiān)測。

4）路基不均勻沉降促使道床產(chǎn)生應(yīng)力集中現(xiàn)象，道床應(yīng)力增大；隨著沉降波長和幅值的增加，進一步導(dǎo)致道床應(yīng)力集中范圍外延，長期運營將導(dǎo)致道砟破碎，道床板結(jié)；期望為工務(wù)維修“科學(xué)修，精確修”提供理論依據(jù)。