鐵路碎石道砟靜態壓碎行為數值模擬

張 徐， 趙春發， 翟婉明

(西南交通大學牽引動力國家重點實驗室，四川 成都610031)

大量現場觀測表明，在列車荷載長期往復作用 下，道床內部的碎石道砟會出現不同程度的破碎與粉化現象［1］. 道砟破碎和粉化會導致道床臟污與板結，使道床彈性下降，并增強輪軌動力作用，加快道床宏觀殘余變形和軌道結構沉降，進而惡化軌面幾何狀態，加劇輪軌動力作用.如此惡性循環，可能導致軌道結構狀態超過安全使用限度，引發行車安全事故.我國既有鐵路經過6 次大提速后，干線旅客列車最高時速達到200 km/h 以上;2008 年以來，建成并開通時速200 km/h 以上高速鐵路1 萬余km(截止2013 年底).伴隨我國鐵路提速與高速化發展，道砟破碎、粉化導致的有砟軌道道床病害問題日益突出，由道床病害引起的軌道養護維修工作量顯著增加，這需要深入研究碎石道砟及散體道床的力學行為與劣化機理，并在掌握道砟破壞和道床變形基本規律的基礎上，發展、完善提速及高速有砟軌道結構設計、養護維修理論與方法.

道砟破碎與道砟顆粒強度及其承受的荷載有關，也與道砟材質及其內部應力分布、微裂縫缺陷有關［2］，研究單個道砟破碎的細觀力學行為有利于揭示道砟顆粒間相互接觸導致破碎的機理.國內外學者通過靜態壓碎試驗研究了道砟顆粒的強度特征，如McDowell、蘇勇等分別測試了粒徑10 ～50 mm 道砟的靜態壓碎強度，指出單顆粒道砟靜態壓碎的特征強度服從Weibull 分布［3-5］.

室內試驗可以獲得道砟的破碎強度，結果也可用于數值模型的驗證分析，但現有測試技術不能獲得加載過程中道砟的內部應力，因此，通過室內試驗研究尚不能完全揭示道砟的破碎機理.

用離散元法模擬道砟的靜態壓碎行為，可以獲得壓碎過程中道砟內部的應力分布與裂紋擴展進程，便于研究道砟顆粒破碎及道床累積變形機理.Guerrero 等用PFC2D 軟件建立了二維軌枕-道床模型，研究了循環荷載作用下道砟破碎對軌道永久變形的影響，發現即使只有少量顆粒發生破碎，也會嚴重增大軌道永久變形［6］;Lim 等采用黏結球形單元，根據加載時單元間黏結鍵失效來模擬單顆粒道砟的靜態壓碎行為［7］;Ergenzinger 等采用球形單元填充多面體的方法模擬單顆粒道砟的壓碎行為［8］.

道砟顆粒靜態壓碎行為的上述離散元模擬研究中，均簡化了道砟顆粒的形狀，也沒有關注壓碎過程中道砟內部應力分布與裂紋擴展的演化.

為深入了解碎石道砟靜態壓碎強度及其壓碎過程中的力學行為，我們使用三維激光掃描儀，獲得了一些高速鐵路特級道砟的真實形態，采用PFC3D 軟件構建道砟靜態壓碎的離散元模型，對碎石道砟靜態壓碎行為進行數值模擬，重點分析了道砟的壓碎強度、壓碎過程中內部應力分布、演化以及裂紋擴展過程.

1 道砟靜態壓碎的離散元分析模型

1.1 道砟幾何形態模擬

碎石道砟具有形狀不規則、棱角分明、表面紋理粗糙的形態特征. 通常可以采用數碼相機、攝影機或CT 設備等提取道砟的幾何形態［9-11］，但一般僅能獲得道砟的二維或部分三維形態特征.為準確獲得道砟真實的三維形態，我們使用三維激光掃描儀(圖1)，獲得了道砟表面點云.圖2 為單顆道砟的實物照片和掃描獲得的三維表面點云，由621 035 個表面點構成.

圖1 單個道砟的三維激光掃描Fig.1 3-D laser scanning of a ballast

圖2 道砟實物和道砟表面點云Fig.2 A ballast and its surface point cloud

利用掃描獲得的道砟三維幾何圖形，采用成密排六方結構的球形單元填充表面點云圍成的區域構建道砟的離散元模型.圖3 為分別采用半徑r =1，2，5 mm 的球形單元構建的道砟的離散元模型.由圖3 可見，球形單元半徑越小，道砟離散元模型的幾何形態越真實.顯然，對于同一道砟，采用的球形單元半徑越小，模擬精度越高，但所需單元數量也越多，這將導致模型的求解規模越大，計算時間越長.因此，為了兼顧模擬精度與計算效率，經試算對比后，采用半徑為2 mm 的球形單元.

圖3 道砟離散元模型Fig.3 Discrete element models for a ballast

1.2 接觸模型與黏結斷裂準則

道砟材料的本構行為通過球形單元間的接觸與黏結實現. 采用線彈性接觸本構模型［12］計算單元間的接觸力，加載前在相互接觸的球形單元間定義接觸黏結，黏結斷裂準則為:

式中:Fn和Fs分別為法向和切向接觸力標量，法向接觸力以受壓為正;φn和φs分別為法向和切向接觸黏接強度.

當單元間接觸力均不滿足式(1)和(2)時，黏結不會發生斷裂;當單元間接觸力滿足式(1)或(2)中之一時，接觸黏結將發生斷裂，法向和切向接觸黏接強度均失效，此后兩單元若再發生接觸，其法向和切向接觸力不再受式(1)和(2)限制. 加載過程中，模型中黏結斷裂的數量逐漸增加，可以直觀模擬道砟在載荷作用下內部裂紋萌生、擴展直至破碎的過程.

1.3 模型參數選取

為了使離散元模型能獲得可靠的宏觀響應，一般需要根據材料的宏觀特性來標定單元顆粒的細觀參數［13-15］.將兩相互接觸的球形單元等效為長度為兩接觸球體半徑之和的彈性梁，根據拉壓和抗彎剛度等效的原則，可以建立材料宏觀特性與顆粒細觀參數的關系［12］.對于花崗巖材料的道砟，參考試驗獲得的花崗巖的宏觀特性參數［16-17］，計算出離散元模型的微觀參數值，見表1.

表1 離散元模型的微觀參數Tab.1 Micro parameters in the discrete element models

1.4 加載方式與加載速率

以黏結的球形單元填充掃描獲得的碎石道砟的三維幾何形狀，以相互平行的兩剛性平板模擬試驗機壓頭，建立碎石道砟靜態壓碎的離散元模型，見圖4.加載時下加載板固定，上加載板以恒定的速率v 向下加載板移動，同時記錄施加在上加載板的載荷大小和位移量. 通常，試驗機壓頭的剛度遠高于巖石剛度，而且金屬壓頭較為光滑，故取加載平板的剛度為5 GN/m，摩擦因數為0.1. 此外，取模型中局部慣性阻尼系數為0.7，忽略材料的粘性阻尼.

圖4 道砟靜態壓碎的離散元模型Fig.4 Discrete element model for static crushed behavior of a ballast

對同一道砟，分別以v =0.1，1.0，10.0 mm/s的加載速率模擬其靜態壓碎過程.數值模擬結果表明，3 種加載速率下的載荷-位移歷程基本重合，載荷峰值差異極小;加載速率為0.1 和1.0 mm/s 時，道砟破碎過程中的單元黏結斷裂數N 十分接近;但是，當加載速率提高到10.0 mm/s 時，單元黏結斷裂數明顯降低. 可見，過大的加載速率會引起一定的慣性沖擊，導致道砟壓碎過程呈現出不同的破壞方式，故下文計算中采用v =1.0 mm/s 的加載速率.

2 碎石道砟靜態壓碎強度統計分析

McDowell 等根據加載時施加在道砟上的載荷F 和上、下加載平板之間的距離d 定義道砟的特征應力σ 為［3］:

當F 取加載過程中的載荷峰值，d 取道砟破壞時加載平板間的距離時，將由式(3)得到的特征應力定義為碎石道砟靜態壓碎特征強度.

與巖石強度的點載荷試驗結果一樣，道砟靜態壓碎試驗獲得的特征強度也呈現很大的離散性.McDowell、蘇勇等的研究表明，道砟靜態壓碎特征強度服從Weibull 概率分布［3-5］，即特征應力為σ時道砟的殘存概率Ps滿足:

式中:σ0是殘存概率為37%時的特征應力;m 為Weibull 系數.

根據高速鐵路特級碎石道砟的級配要求，采用離散元法，模擬分析了粒徑分別為30 ～40、40 ～50和50 ～60 mm 共27 顆道砟的靜態壓碎行為.經統計分析得到不同粒徑道砟特征強度的殘存概率分布，見圖5(x 表示橫坐標，y 表示縱坐標，R 為相關系數). 圖5 中還給出了相應的線性擬合公式、擬合曲線以及殘存概率為37%時的Weibull統計特征強度.

圖5 不同粒徑道砟靜態壓碎特征強度的Weibull 分布Fig.5 Weibull survival probability of railway ballast with different particle sizes

ln(ln(1/Ps))與ln σ 的高度線性相關性表明，用Weibull 分布能夠較好地描述碎石道砟靜態壓碎時特征強度的分布規律.

從圖5 還可見:

(1)粒徑30 ～40、40 ～50 和50 ～60 mm 的碎石道砟的Weibull 系數m 分別為2.348 8、3.336 4和1.279 4，對應的殘存概率為37%的特征強度分別為37.48、25.39 和20.23 MPa，與試驗結果［3-4］接近，說明建立的離散元模型合理，計算結果有較高的可信度.

(2)粒徑越大，道砟殘存概率為37%的特征強度越小.

進一步分析表明，碎石道砟靜態壓碎殘存概率為37%的特征強度與其平均粒徑具有近似的線性相關性.

3 碎石道砟靜態壓碎的力學行為

碎石道砟靜態壓碎時，不僅壓碎強度離散性大，壓碎過程中道砟內部應力分布與裂隙擴展過程也存在差異.但是，道砟從初始加載至最終整體破壞，其靜態壓碎過程基本包含4 個階段:顆粒轉動、局部破碎、彈性響應和整體破壞階段［18］.需要指出的是，在某個確定道砟的靜態壓碎過程中，有可能多個階段不出現，也有可能某個階段在加載過程中多次出現，這與道砟的形態特征、加載壓頭與道砟接觸位置等有關.

建立某確定道砟的離散元模型，模型包含658 個半徑為2 mm 的球形單元，所有單元間初始接觸黏結數為3 162 個.為揭示該道砟靜態受壓時的破碎機理，從載荷-位移響應、內部接觸力鏈演化與裂隙擴展三方面分析道砟靜態壓碎過程中不同階段的力學行為.

圖6 給出了該道砟加載前和壓碎后的離散元模型.圖6(b)表明，該道砟經靜態壓縮后發生整體破壞，分裂成2 個主要碎塊.

圖6 道砟在靜態壓碎前、后的離散元模型Fig.6 Discrete element models of a ballast when it is intact and crushed

3.1 載荷-位移響應

圖7 為道砟從初始加載到整體破壞的載荷-位移(F-s)響應曲線.初始加載時，由于道砟與下加載板的接觸表面不平整，載荷作用下道砟發生轉動，此時載荷處于極低水平;道砟轉動至穩定位置后，與上、下加載板的作用點逐漸固定，處于較穩定的受力狀態，道砟變形量線性增大;當加載位移增大到0.396 mm 時，道砟與加載板接觸點附近局部產生微小破碎，道砟發生第2 次輕微轉動，此時載荷急劇減小至接近于0，表明上加載板與道砟基本上處于脫離接觸的狀態;當加載位移以恒定速率繼續增大時，上加載板與道砟重新接觸，道砟再次進入穩定受力狀態，載荷隨位移增大線性增大;最后，隨載荷和位移的進一步增大，道砟內部部分應力超過黏結強度，內部產生裂紋并迅速擴展，導致道砟整體破壞，載荷瞬間減小至0.

圖7 道砟靜態壓碎過程中載荷-位移響應歷程Fig.7 Load-displacement history in ballast crush process

3.2 道砟內部力鏈分布及其演變

為掌握靜態壓碎過程中道砟內部發生的細觀力學行為，圖8 給出了不同階段道砟內部接觸力鏈的分布及其演變過程. 圖8 中，連接兩接觸單元球心的實線表示單元間的接觸力，紅色表示壓力，綠色表示拉力，線條越粗則接觸力幅值越大.

圖8(a)為初始穩定階段道砟內部接觸力鏈的分布，表明道砟與兩加載板間形成了幾個穩定的接觸點，這些點處的接觸力明顯高于道砟內部單元的接觸力，符合道砟受壓時載荷從表面若干個接觸點向內部傳遞的自然規律.

圖8(b)為載荷達到6.23 kN 時，道砟內部接觸力鏈的分布. 結合載荷-位移響應曲線(圖7)可知，當載荷從0 增大到6.23 kN 時，道砟與下加載板之間某個接觸點附近發生局部破碎，打破了道砟與加載板之間的穩定接觸狀態，道砟發生第2 次轉動.

道砟第2 次轉動后，加載位移繼續增大，加載板上的載荷從0 開始重新增大.圖8(c)給出了載荷重新達到6.14 kN 時道砟內部接觸力鏈的分布.對比圖8(b)和8(c)可見，道砟與下加載板的接觸點位置明顯變化，且道砟內部力鏈的分布也發生了較大變化.

圖8 道砟靜態壓碎過程中內部接觸力鏈的演化Fig.8 Evolution of force chain distribution in the crush process

圖8(d)～8(f)為道砟從承受最大荷載到瞬間破碎時內部力鏈的分布.圖8(d)顯示，當載荷達到23.77 kN 時，道砟下部局部區域集中出現單元接觸拉力，相當數量的單元接觸力超過其黏接強度，道砟內部出現了一定規模的黏結斷裂，形成微裂隙并迅速沿主力鏈方向擴展，最終導致道砟整體破壞.圖8(f)表明，當道砟破碎為2 個碎塊后，上加載板與道砟脫離接觸，2 個碎塊之間不再有力鏈連接，碎塊內部的力鏈分布基本恢復到自然狀態.

3.3 道砟內部黏結斷裂及其分布

圖9 為靜態壓碎過程中道砟內部單元接觸黏結斷裂數隨載荷的變化曲線.根據道砟內部單元接觸黏結斷裂數隨載荷增大的快慢，可將道砟從內部裂隙萌生到整體破壞過程分為3 個階段:裂隙萌生與聚集、裂隙急劇擴展和道砟整體破壞. 在裂隙萌生與聚集階段，當載荷從0 增大到23.77 kN 時，雖然道砟內部黏結斷裂數由0 增加到40 個，道砟內部出現若干微裂紋，但總體上載荷與加載位移近似成正比關系(見圖7)，表明道砟整體上處于彈性變形，并具有相當大的承載能力;在裂隙急劇擴展階段，內部黏結斷裂數瞬間增加到168 個，道砟內部已經出現裂縫，同時伴隨著載荷減小，表明道砟承載能力開始下降;在道砟整體破壞階段，黏結斷裂數快速增加到267 個，載荷急劇下降，道砟已經失去承載能力.

圖9 靜態壓碎過程中內部黏結斷裂數隨載荷的變化Fig.9 Number of broken contact bonds vs. load in the crush process

圖10 為靜態壓碎過程中不同階段道砟內部黏結斷裂的分布(紅色方塊表示黏結斷裂的位置).圖10(a)顯示，當載荷首次達到6.23 kN 時，由于道砟與加載板接觸點處存在應力集中，下加載板接觸點附近出現1 處黏結斷裂. 當道砟發生輕微轉動、重新穩定受載后，載荷隨加載位移持續增大，道砟內部黏結斷裂并沒有迅速增加，主要以整體彈性變形來承受外部載荷.直到載荷達到13.60 kN 時，道砟內部出現第2 處黏結斷裂，位于上加載板與道砟接觸點附近(圖10(b)). 當載荷進一步增大到18.86 kN 時(圖10(c))，道砟與上、下加載板接觸點附近黏結斷裂數增加到8 個，這一區域的集中斷裂使得道砟出現較大的局部壓縮變形，道砟與加載板接觸狀態突然發生輕微變化，載荷在較短時間內下降到16.64 kN，且應力集中區域黏結斷裂出現連鎖反應，黏結斷裂數快速增加到15 個.

圖10(d)顯示，當載荷達到峰值23.77 kN 時，道砟內部黏結斷裂數為40 個，且大多集中在沿加載方向的某一平面或其附近區域.該區域大致位于道砟與下加載板接觸面中部的上方，承受剪切力和較大的彎矩.此時，道砟在該區域開始形成集中的微裂隙，并沿加載方向迅速擴展，當加載位移進一步增大時載荷開始下降. 圖10(e)表明，當載荷由峰值降到22.26 kN 時，黏結斷裂數由40 個急劇增加到168 個，說明道砟在整體破壞前期裂隙擴展十分迅速.比較圖10(d)和10(e)可知，大量新增黏結斷裂多出現在黏結斷裂面的擴展面上，且主要沿加載方向向上迅速發展，形成明顯的斷裂面. 這一斷裂面在圖10(e)中更為清晰，而且可以發現，在斷裂面處單元間已經不再接觸.

圖10 道砟靜態壓碎過程中黏結斷裂分布的演化Fig.10 Evolution of the distribution of broken bonds in the crush process

繼續加載時，道砟整體破壞.圖10(f)顯示，在較小的載荷作用下，道砟內部黏結斷裂數快速增加到267 個，黏結斷裂面已沿加載方向貫穿整個道砟，使道砟迅速劈裂成2 個碎塊.

4 結 論

本文中采用激光掃描儀獲得鐵路碎石道砟的真實三維幾何形態，建立了道砟靜態壓碎的離散元模型，對道砟靜態壓碎行為進行了數值模擬，通過統計道砟靜態壓碎的特征強度，分析了道砟壓碎過程中載荷-位移響應、內部力鏈分布與演化以及單元接觸黏結斷裂的發展過程，得到以下結論:

(1)道砟靜態壓碎的特征強度服從Weibull分布，道砟平均粒徑增大時，其殘存概率為37%的特征強度會降低，這與已有試驗結果一致. 說明所建立的道砟靜態壓碎離散元模型正確，模型參數選取合理，可供鐵路散體道床細、宏觀力學行為及狀態劣化模擬研究參考.

(2)道砟尖銳棱角和表面不平整容易引起道砟發生翻轉和局部壓碎，導致道砟外部接觸狀態發生微弱變化，載荷也出現短暫回落;道砟穩定受載時，載荷與加載位移近似成正比;道砟內部裂紋快速擴展和整體破壞階段，載荷急劇減小.

(3)道砟與加載板接觸點處的接觸力明顯高于道砟內部單元的接觸力;道砟外部接觸狀態發生變化會導致其內部接觸力鏈分布明顯變化;道砟穩定受載時，內部接觸力鏈分布穩定，隨加載位移增大，內部黏結斷裂數逐漸增大;當載荷達到峰值后，內部黏結斷裂數急劇增加，道砟劈裂破碎.

致謝:西南交通大學科技創新項目(2682014CX043).參考文獻:

［1］ 曾樹谷. 鐵路散粒體道床［M］. 北京:中國鐵道出版社，1997:30-57.

［2］ AURSUDKIJ B. A laboratory study of railway ballast behaviour under traffic loading and tamping maintenance［D］. Nottingham:University of Nottingham，2007.

［3］ MCDOWELL G R，LIM W L，COLLOP A. Measuring the strength of railway ballast［J］. Ground Engineering，2003，36(1):25-28.

［4］ LIM W L，MCDOWELL G R，COLLOP A C. The application of Weibull statistics to the strength of railway ballast［J］. Granular Matter，2004，6(4):229-237.

［5］ 蘇勇. 含沙鐵路道碴力學行為研究［D］. 大連:大連理工大學工程力學系，2011.

［6］ LOBO G S，VALLEJO L E. Discrete element method analysis of railtrack ballast degradation during cyclic loading［J］. Granular Matter，2006，8(3/4):195-204.

［7］ LIM W L， MCDOWELL G R. Discrete element modelling of railway ballast［J］. Granular Matter，2005，7(1):19-29.

［8］ ERGENZINGER C，SEIFRIED R，EBERHARD P. A discrete element model predicting the strength of ballast stones［J］. Computers ＆ Structures，2012，108:3-13.

［9］ FERELLEC J F，MCDOWELL G R. A method to model realistic particle shape and inertia in DEM［J］.Granular Matter，2010，12(5):459-467.

［10］ 杜欣，曾亞武，高睿. 基于CT 掃描的不規則外形顆粒三維離散元建模［J］. 上海交通大學學報，2011，45(5):711-715.DU Xin，ZENG Yawu，GAO Rui. 3D modelling of irregular shape particles for discrete element method based on X-ray tomography［J］. Journal of Shanghai Jiaotong University，2011，45(5):711-715.

［11］ ANOCHIE B J K，KOMBA J，TUTUMLUER E.Aggregate surface areas quantified through laser measurements for South African asphalt mixtures［J］.Journal of Transportation Engineering，2012，138(8):1006-1015.

［12］ Itasca Consulting Group，Inc. Particle flow code in 3 dimensions (PFC3D)［M］. 4th ed. Minneapolis，MN:Itasca Consulting Group，Inc.，2008:3-11-3-21.

［13］ STAHL M， KONIETZKY H. Discrete element simulation of ballast and gravel under special consideration of grain-shape，grain-size and relative density［J］. Granular Matter，2011，13(4):417-428.

［14］ WANG Y，TONON F. Modeling Lac du Bonnet granite using a discrete element model［J］. International Journal of Rock Mechanics and Mining Sciences，2009，46(7):1124-1135.

［15］ WANG Y，TONON F. Calibration of a discrete element model for intact rock up to its peak strength［J］. International Journal for Numerical and Analytical Methods in Geomechanics，2010，34(5):447-469.

［16］ POTYONDY D O，CUNDALL P A. A bonded-particle model for rock［J］. International Journal of Rock Mechanics and Mining Sciences，2004，41(8):1329-1364.

［17］ 王渭明，楊更社，張向東，等. 巖石力學［M］. 徐州:中國礦業大學出版社，2010:36-48.

［18］ CAVARRETTA I，O'SULLIVAN C. The mechanics of rigid irregular particles subject to uniaxial compression［J］. Geotechnique，2012，62(8):681-692.