張 徐， 趙春發(fā)， 翟婉明

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

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

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

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

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

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

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

1 道砟靜態(tài)壓碎的離散元分析模型

1.1 道砟幾何形態(tài)模擬

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

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

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

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

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

1.2 接觸模型與黏結(jié)斷裂準則

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

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

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

1.3 模型參數(shù)選取

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

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

1.4 加載方式與加載速率

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

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

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

2 碎石道砟靜態(tài)壓碎強度統(tǒng)計分析

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

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

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

式中:σ0是殘存概率為37%時的特征應(yīng)力;m 為Weibull 系數(shù).

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

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

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

從圖5 還可見:

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

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

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

3 碎石道砟靜態(tài)壓碎的力學行為

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

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

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

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

3.1 載荷-位移響應(yīng)

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

圖7 道砟靜態(tài)壓碎過程中載荷-位移響應(yīng)歷程Fig.7 Load-displacement history in ballast crush process

3.2 道砟內(nèi)部力鏈分布及其演變

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

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

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

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

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

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

3.3 道砟內(nèi)部黏結(jié)斷裂及其分布

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

圖9 靜態(tài)壓碎過程中內(nèi)部黏結(jié)斷裂數(shù)隨載荷的變化Fig.9 Number of broken contact bonds vs. load in the crush process

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

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

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

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

4 結(jié) 論

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

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

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

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

致謝:西南交通大學科技創(chuàng)新項目(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 掃描的不規(guī)則外形顆粒三維離散元建模［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］. 徐州:中國礦業(yè)大學出版社，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.