孫珊珊， 嚴(yán) 穎， 趙春發(fā)， 季順迎

(1. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023； 2. 大連交通大學(xué) 土木與安全工程學(xué)院，遼寧 大連 116028； 3. 西南交通大學(xué) 牽引動(dòng)力國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610031)

有砟軌道是當(dāng)前鐵路軌道的主要結(jié)構(gòu)形式，在鐵路運(yùn)輸向高速和重載方向發(fā)展過程中，其力學(xué)性能不斷改進(jìn)[1-2]。鐵路道床在列車往復(fù)荷載作用下因道砟顆粒的重新排列和破碎粉化發(fā)生沉降，導(dǎo)致整體彈性的降低，影響其使用性能[3]。

相關(guān)學(xué)者對材料的力學(xué)性能開展了系統(tǒng)研究，主要集中在累積沉降量和有效彈性等方面。道砟材料的三軸試驗(yàn)結(jié)果表明，在密實(shí)過程中道砟材料減震和排水性能均有所降低[4-5]。在高頻重載條件下，道砟材料的沉降和道砟顆粒破壞現(xiàn)象尤為突出[6]。在往復(fù)荷載作用下，道砟顆粒的破碎是引起道床變形的重要原因[7]。將鐵路道床假設(shè)為連續(xù)體，用有限元方法可分析其動(dòng)力特性和變形規(guī)律[8-9]，但僅限小變形行為，且不能描述道砟顆粒的細(xì)觀受力特征。顆粒的幾何形狀及級配對其宏觀力學(xué)行為的影響顯著[10-11]。

為真實(shí)模擬道砟材料的幾何形態(tài)和力學(xué)行為，基于非連續(xù)材料的離散單元模型迅速發(fā)展[12-14]。在顆粒材料的離散元數(shù)值模擬中，采用最普遍的是球形顆粒單元。對于非規(guī)則形態(tài)的道砟顆粒，球形顆粒很難有效地模擬道砟材料的力學(xué)行為[15]，可采用球形顆粒的不同黏接或鑲嵌組合方式構(gòu)造相應(yīng)的道砟顆粒[14,16-17]，模擬道砟顆粒在外荷載作用下的破碎現(xiàn)象，對有砟道床的沉降特性進(jìn)行離散元分析[7,17-18]。在球形顆粒的組合中，通過增加球形顆粒單元的數(shù)量可更精確地模擬非規(guī)則顆粒，甚至可采用上千個(gè)球形顆粒單元構(gòu)造一個(gè)道砟塊石[16]。但顆粒數(shù)目的增加急劇降低離散元的計(jì)算效率，不利于開展大規(guī)模計(jì)算。

為模擬非球形顆粒的幾何形態(tài)，發(fā)展了不同的離散單元模型。采用二次或超二次曲面方程可定義不同形態(tài)的顆粒單元[19-20]。對類似道砟材料的多棱邊顆粒，多面體單元能夠合理描述顆粒的幾何形態(tài)以及顆粒間的相互作用，準(zhǔn)確模擬其宏觀力學(xué)特征[21]。多面體單元的作用力計(jì)算需要繁瑣的接觸判斷和鄰居搜索，嚴(yán)重影響其計(jì)算效率和工程應(yīng)用性。采用Minkowski Sum方法構(gòu)造擴(kuò)展多面體顆粒模型，準(zhǔn)確快速地模擬非規(guī)則顆粒間的相互作用。在離散單元法的數(shù)值研究中，這種非規(guī)則顆粒模型已成功用于二維的擴(kuò)展多邊形[22-23]和三維的擴(kuò)展多面體計(jì)算[24-26]。擴(kuò)展多面體單元的主要特點(diǎn)是在多面體單元的棱邊及外表面設(shè)置若干個(gè)球形顆粒，通過改變球形顆粒的直徑靈活控制多面體棱邊及棱角的尖銳度[27]。該模型可有效簡化多面體單元的接觸判斷過程，在滿足計(jì)算精度的條件下有效地提高離散元計(jì)算效率。

本文采用擴(kuò)展多面體模型對有砟道床在往復(fù)荷載作用下的動(dòng)力過程進(jìn)行離散元分析，確定其累積沉降量和有效模量隨荷載往復(fù)次數(shù)的演化規(guī)律。

1 道砟顆粒的擴(kuò)展多面體離散元模型

基于Minkowski Sum 方法，道砟顆粒的擴(kuò)展多面體單元由基本多面體單元和擴(kuò)展球體單元疊加而成。Minkowski Sum方法是處理復(fù)雜圖像和構(gòu)造復(fù)雜顆粒模型的一種有效方法[25-26]。將2個(gè)空間體A、B疊加，即在歐幾里得幾何空間中將一個(gè)空間點(diǎn)集掃過另一個(gè)點(diǎn)集，得到新的集合，其定義式為

A⊕B={x+y|x∈A，y∈B}

( 1 )

式中:A、B分別為2個(gè)空間體；x、y分別為其對應(yīng)的空間坐標(biāo)。

擴(kuò)展單元為球體時(shí)，它將空間上的一個(gè)點(diǎn)擴(kuò)展為球體，將一條線段擴(kuò)展為球形圓柱體，將任一平面擴(kuò)展為由圓柱體包絡(luò)的具有一定厚度的板，將任一多面體擴(kuò)展為具有光滑棱邊及棱角的多面體[28-29]。

本文構(gòu)造的擴(kuò)展多面體單元是將多面體作為基本單元，將球體作為擴(kuò)展單元，將2個(gè)空間體疊加，形成具有光滑棱邊和角點(diǎn)的非規(guī)則形態(tài)單元體，見圖1。

在離散元數(shù)值模擬中，通過改變擴(kuò)展球體的直徑得到不同光滑度的非規(guī)則顆粒單元。設(shè)基本立方體單元和四面體單元邊長均為a，球體顆粒半徑r分別為a/10、a/5時(shí)構(gòu)造的擴(kuò)展多面體單元見圖2。從圖2可以看出球體顆粒半徑對擴(kuò)展多面體單元棱邊和棱角尖銳度的影響。

針對道砟顆粒的非規(guī)則幾何形態(tài)，本文采用擴(kuò)展多面體模型構(gòu)造不同形態(tài)的顆粒單元，見圖3。由于擴(kuò)展多面體單元具有非規(guī)則性，用微元?jiǎng)澐址ù_定其體積，繼而確定其質(zhì)量、形心及對形心的轉(zhuǎn)動(dòng)慣量。在顆粒單元的作用力和運(yùn)動(dòng)計(jì)算中，采用全局坐標(biāo)和局部坐標(biāo)2種坐標(biāo)。全局坐標(biāo)固定在數(shù)值模擬的總體計(jì)算局域上，局部坐標(biāo)固定在每個(gè)單元的形心上并隨顆粒的轉(zhuǎn)動(dòng)而改變。采用歐拉角表征每個(gè)單元局部坐標(biāo)的方位，通過四元數(shù)方法計(jì)算局部坐標(biāo)與全局坐標(biāo)之間的轉(zhuǎn)換矩陣，對顆粒間作用力和轉(zhuǎn)動(dòng)分量進(jìn)行坐標(biāo)轉(zhuǎn)化[30]。

2 鐵路道砟單元接觸判斷

擴(kuò)展多面體道砟顆粒主要由平面、棱邊和棱角3部分組成，它們在道砟材料動(dòng)力過程中發(fā)生相互接觸，主要有3類接觸形式：(1) 棱角-平面、棱角-棱邊及棱角-棱角接觸；(2) 棱邊-棱邊、棱邊-平面接觸；(3)平面-平面接觸。本文選取有代表性的棱角-平面、棱邊-棱邊及平面-平面接觸進(jìn)行接觸判斷。

2.1 棱角-平面接觸

道砟材料的棱角與平面間的作用狀態(tài)見圖4(a)。擴(kuò)展多面體單元中道砟棱角為球形顆粒，其與平面的接觸判斷轉(zhuǎn)化為球體與平面的接觸判斷，見圖4(b)。已知球體球心O1和平面中心O2的位置矢量，由幾何關(guān)系可確定球體與平面點(diǎn)的重疊量。

確定球心O1到平面的距離矢量

Δ=n(O21·n)

( 2 )

式中：O21為O2到O1的位置矢量；n為平面的法向單位矢量。

O21在平面上的分量

P21=O21-Δ

( 3 )

式中：P21為O2到球心O1在平面上的投影P的位置矢量。

判斷P點(diǎn)是否在平面內(nèi)。P點(diǎn)不在平面內(nèi)，棱角球體與平面不發(fā)生接觸；P點(diǎn)在平面內(nèi)時(shí)，球體到平面的最短距離

δ=|Δ|-2R

( 4 )

式中：R為擴(kuò)展多面體單元選用的球體顆粒半徑；δ為發(fā)生接觸時(shí)2個(gè)顆粒單元間的法向重疊量。

δ<0時(shí)發(fā)生接觸，否則不接觸。發(fā)生接觸時(shí)，接觸力法線與Δ平行，接觸點(diǎn)分別為Pc1、Pc2。對棱角-平面的接觸位移進(jìn)行計(jì)算，確定接觸作用力。

2.2 棱邊-棱邊接觸

2.3 平面-平面接觸

擴(kuò)展多面體單元中平面-平面的接觸判斷見圖6。為簡化計(jì)算，將平面-平面接觸轉(zhuǎn)化為接觸面上的多點(diǎn)接觸。圖6中，將2個(gè)平面間的接觸力通過4個(gè)點(diǎn)接觸進(jìn)行計(jì)算，其中接觸點(diǎn)A、C為棱角-平面接觸，B、D為棱邊-棱邊接觸。4個(gè)接觸點(diǎn)的接觸力確定后，通過各接觸力對顆粒質(zhì)心的力矩確定接觸面間的力矩。該方法具有模型簡單、計(jì)算效率高的特點(diǎn)。

3 道砟單元間的接觸模型

擴(kuò)展多面體單元間接觸模式不同，接觸面積差異大，導(dǎo)致單元間的接觸剛度和阻尼系數(shù)不同。計(jì)算擴(kuò)展多面體單元的接觸力時(shí)需考慮接觸點(diǎn)處單元的曲率半徑和接觸面積。

3.1 道砟單元間的接觸力模型

擴(kuò)展多面體道砟單元間的接觸力主要取決于2個(gè)接觸單元間的重疊量和相對速度。在不同接觸模式下，接觸單元間的法、切向力分別為

( 5 )

( 6 )

式( 5 )、式( 6 )中各參數(shù)與擴(kuò)展多面體單元的接觸模式密切相關(guān)。

在接觸的切線方向，依據(jù)Mohr-coulomb摩擦定律，有

( 7 )

式中:μ為道砟顆粒材料的摩擦系數(shù)。

擴(kuò)展多面體單元的棱角-棱角接觸中，接觸方式為擴(kuò)展球體的接觸。采用Hertz非線性接觸模型計(jì)算其接觸力[31-32]，其法向接觸剛度為

( 8 )

式中：G、ν分別為道砟材料的剪切模量和泊松比；R*為接觸球體的有效半徑，R*=R1R2/(R1+R2)，其中R1、R2為2個(gè)擴(kuò)展球體的曲率半徑。

切向剛度與法向剛度成線性關(guān)系，Kt=βKn。系數(shù)β由道砟材料的泊松比確定，β=3(1-ν)/(2-ν)。

計(jì)算接觸單元間法向和切向黏滯力時(shí)，法向、切向黏滯力系數(shù)分別為

( 9 )

(10)

3.2 時(shí)間步長的確定

對于2個(gè)球體單元，計(jì)算時(shí)間步長采用簡諧振動(dòng)模型確定[33]。忽略法向黏滯力的影響，基于簡諧振動(dòng)方程得最小振動(dòng)周期

(11)

式中:Mmin為最小顆粒質(zhì)量。

按2個(gè)球體單元的接觸模式確定時(shí)間步長

Δt=Tmin/40

(12)

4 道砟材料動(dòng)力特性的擴(kuò)展多面體離散元分析

4.1 道砟箱試驗(yàn)的離散元模型

采用道砟箱試驗(yàn)對道砟材料在往復(fù)荷載下的動(dòng)力特性進(jìn)行離散元分析。道砟顆粒材料的隨機(jī)生成過程中，調(diào)整道砟顆粒的長寬比，生成多種顆粒單元。道砟顆粒最大粒徑在36～63 mm均勻分布，所有擴(kuò)展多面體道砟顆粒的擴(kuò)展球體半徑相同。主要計(jì)算參數(shù)見表1。

表1 道砟動(dòng)力特性離散元模擬中的主要計(jì)算參數(shù)

考慮有砟道床的真實(shí)結(jié)構(gòu)和軌枕距離，參考文獻(xiàn)[15]中采用球體組合顆粒對道砟材料動(dòng)力特性的離散元分析，取計(jì)算域長L=700 mm，寬B=300 mm；軌枕長Ls=300 mm，寬Bs=300 mm，質(zhì)量Ms=34 kg。為模擬有砟道床在無限長度的動(dòng)力特性，在長度方向采用周期邊界條件。有砟道床結(jié)構(gòu)的道砟厚度一般不小于350 mm，取Hs=370 mm。

將道砟顆粒自由下落到道砟箱內(nèi)，在蓋板上施加頻率3 Hz、大小0～3 kN的循環(huán)荷載，持續(xù)3 s以壓實(shí)。將軌枕放到道床上，對軌枕以低速率加載至3 kN，完成道砟顆粒的初始排列。在軌枕上施加正弦往復(fù)荷載，最大、最小值分別為40、3 kN，荷載頻率3 Hz。只考慮軌枕在外荷載和道砟反力的作用下的豎直位移，忽略其水平位移和扭轉(zhuǎn)運(yùn)動(dòng)。

4.2 往復(fù)荷載下道砟材料的沉降量及變形模量

對道砟材料在往復(fù)荷載作用下的動(dòng)力過程進(jìn)行數(shù)值分析，獲得90 s內(nèi)軌枕位移和相應(yīng)的往復(fù)荷載隨時(shí)間的變化，見圖8。從軌枕位移u的變化趨勢可以發(fā)現(xiàn)，軌枕的豎直沉降量隨往復(fù)荷載呈周期性變化，其累計(jì)沉降量隨加載次數(shù)的增加不斷增大，增長趨勢趨于平緩。對78～80 s的計(jì)算結(jié)果進(jìn)行放大，可清晰地看出軌枕位移和荷載的對應(yīng)關(guān)系，見圖8(b)。軌枕的上下震動(dòng)與往復(fù)荷載在頻率和相位上均保持一致。

在90 s內(nèi)的270個(gè)加載周期內(nèi)，外加荷載與軌枕沉降量的對應(yīng)關(guān)系見圖9。從圖9可以看出，最初幾次加載時(shí)，軌枕的位移明顯較大，因?yàn)榈理牟牧系某跏寂帕邢鄬κ杷?，在外荷載的作用力下達(dá)到密實(shí)排列。隨著加載次數(shù)的增加，道砟不斷密實(shí)，位移變形曲線越來越密集。在往復(fù)荷載作用下，道砟材料的累計(jì)沉降量趨于穩(wěn)定。在一個(gè)加載周期內(nèi)，加載過程中位移的增長速率比卸載時(shí)的位移降低速率快。圖9體現(xiàn)了道砟材料變形的2個(gè)特性：(1)彈性變形，即每次加卸載過程中軌枕的豎向位移變化量，體現(xiàn)有砟道床的彈性性能，主要由道砟顆粒的彈性變形引起；(2)塑性變形，即在多次往復(fù)加卸載過程中軌枕的累計(jì)沉降量，主要由道砟顆粒的重新排列引起。

為分析道床在往復(fù)荷載下彈性性能的演化趨勢，通過每次往復(fù)荷載中軌枕的豎直位移變化量確定道床的有效剛度Ke

Ke=ΔF/Δu

(13)

式中：ΔF為往復(fù)荷載的幅度；Δu為一次加卸載過程中軌枕的豎向位移。

往復(fù)荷載過程中道砟材料有效剛度和軌枕累計(jì)沉降量隨加載次數(shù)的變化見圖10。從圖10(a)可以看出，隨著加載次數(shù)的增加，道床有效剛度不斷增加。在加載的前幾個(gè)周期，由于道床剛度的幾何狀態(tài)尚不穩(wěn)定，其剛度的波動(dòng)性較大。特別是第一個(gè)加載周期，有效剛度較低，為7 kN/mm。隨著加載次數(shù)的增加，道床結(jié)構(gòu)趨于穩(wěn)定，其有效剛度增加到11 kN/mm。圖10(b)中，umin和umax分別為每次往復(fù)荷載達(dá)到最小和最大值時(shí)所對應(yīng)的軌枕沉降量。道床的累計(jì)沉降量有一個(gè)先快速沉降后趨于穩(wěn)定的過程。有效剛度和累計(jì)沉降量在一定程度上反映了有砟道床在往復(fù)荷載作用下的密實(shí)過程，分別表征了道床彈性和永久變形的變化規(guī)律。

5 結(jié)語

為描述鐵路有砟道床中道砟顆粒的非規(guī)則幾何形態(tài)，本文基于Minkowski Sum方法，將多面體模型與球體單元疊加，形成具有非規(guī)則形態(tài)的光滑顆粒單元。將塊體接觸轉(zhuǎn)化為球體單元的接觸，有效簡化了多面體單元間的接觸判斷和接觸力計(jì)算。通過改變擴(kuò)展球體粒徑，構(gòu)造具有不同尖銳度的擴(kuò)展多面體單元，合理模擬非規(guī)則顆粒間相互作用過程。

本文采用擴(kuò)展多面體顆粒單元對道砟材料在往復(fù)荷載下的動(dòng)力過程進(jìn)行了離散元分析，確定了其有效剛度和累計(jì)沉降量隨加載次數(shù)的演化規(guī)律。結(jié)果表明，軌枕位移在隨往復(fù)荷載上下震蕩過程中，有效剛度和累計(jì)沉降量在初始階段波動(dòng)較大，隨道砟材料的密實(shí)度增加而趨于穩(wěn)定。該擴(kuò)展多面體單元能夠合理地模擬道砟材料的動(dòng)力行為，由微觀作用過程揭示相應(yīng)的宏觀演化規(guī)律。在此工作基礎(chǔ)上，今后將對該模型進(jìn)一步改進(jìn)，系統(tǒng)地分析道砟顆粒的級配、摩擦系數(shù)、彈性模量等參數(shù)及道砟在往復(fù)荷載下的破碎和粉化過程對道床力學(xué)性質(zhì)的影響。

參考文獻(xiàn)：

[1] 練松良. 軌道工程[M]. 上海: 同濟(jì)大學(xué)出版社, 2006.

[2] 趙國堂. 高速鐵路無砟軌道結(jié)構(gòu)[M]. 北京: 中國鐵道出版社, 2006.

[3] SUN Q D, INDRARATNA B, NIMBALKAR S. Effect of Cyclic Loading Frequency on the Permanent Deformation and Degradation of Railway Ballast[J]. Geotechnique, 2014，64(9): 746-751.

[4] SUIKER A S J, SELIG E T, FRENKEL R. Static and Cyclic Triaxial Testing of Ballast and Subballast[J].Journal of Geotechnical and Geoenvironmental Engineering, 2005, 131(6): 771-782.

[5] ASCE B I F, THAKUR P K, VINOD J S,et al. Semiempirical Cyclic Densification Model for Ballast Incorporating Particle Breakage[J].International Journal of Geomechanics, 2012, 12(3): 260-271.

[6] LACKENBY J, INDRARATNA B, MCDOWELL G, et al. Effect of Confining Pressure on Ballast Degradation and Deformation under Triaxial Loading[J].Geotechnique,2007,57(6):527-536.

[7] HOSSAIN Z, INDRARATNA B, DARVE F,et al. DEM Analysis of Angular Ballast Breakage under Cyclic Loading[J].Geomechanics and Geoengineering: An International Journal, 2007, 2(3): 175-181.

[8] ZHAI W M, WANG K Y, LIN J H. Modeling and Experiment of Railway Ballast Vibrations[J]. Journal of Sound and Vibration，2004，270(4-5): 673-683.

[9] 董亮, 蔡德鉤, 葉陽升, 等. 列車循環(huán)荷載作用下高速鐵路路基累積變形預(yù)測方法[J].土木工程學(xué)報(bào)，2010, 43(6): 100-108.

DONG Liang, CAI De-gou, YE Yang-sheng, et al. A Method for Predicting the Deformation of High-speed Railway Subgrades under Cyclic Train Loads[J]. China Civil engineering journal, 2010, 43(6): 100-108.

[10] 井國慶, 封坤, 高亮, 等. 循環(huán)荷載作用下道砟破碎老化的離散元仿真[J]. 西南交通大學(xué)學(xué)報(bào), 2012, 47(2): 187-191.

JING Gou-qing, FENG Kun, GAO Liang, et al. DEM Simulation of Ballast Degradation and Breakage under Cyclic Loading[J]. Journal of Southwest Jiaotong University, 2012, 47(2): 187-191.

[11] 肖宏, 高亮, 侯博文. 鐵路道床振動(dòng)特性的三維離散元分析[J]. 鐵道工程學(xué)報(bào), 2009, (9): 14-17.

XIAO Hong, GAO Liang, HOU Bo-wen. Analysis of Ballast Dynamic Behavior with Three-dimensional Discrete Element Method[J]. 2009, (9): 14-17.

[12] MCDOWELL G R, LIM W L, COLLOP A C, et al. Comparison of Ballast Index Tests for Railway Trackbeds[J]. Proceeding of the Institution of Civil Engineers Geotechnical Engineering, 2004, 157(3): 151-161.

[13] 姜衛(wèi)利, 范俊杰. 散粒體道床離散單元法分析[J]. 鐵道學(xué)報(bào), 2001, 23(4): 98-101.

JIANG Wei-li, FAN Jun-jie. Using of Distinct Element Method to Analyze Granular Ballast Bed[J]. Journal of the China Railway Society, 2001,23(4):98-101.

[14] 嚴(yán)穎, 狄少丞, 蘇勇, 等. 風(fēng)沙影響下鐵路道砟變形模量的離散元法數(shù)值分析[J]. 計(jì)算力學(xué)學(xué)報(bào), 2012, 29(3): 339-445.

YAN Ying, DI Shao-cheng, SU Yong, et al.Discrete Element Analysis of Elastic Modulus of Railway Ballasts in Wind with Different Sand Contents[J]. Chinese Journal of Computational Mechanics, 2012, 29(3): 339-445.

[15] LIM W L, MCDOWELL G R. Discrete Element Modeling of Railway Ballast[J]. Granular Matter, 2005, 7(1): 19-29.

[16] LU M, MCDOWELL G R.The Importance of Modeling Ballast Particle Shape in the Discrete Element Method[J]. Granular Matter,2007, 9(1-2): 69-80.

[17] LOBO-GUERRERO S, VALLEJO L. Discrete Element Method Analysis of Rail Track Ballast Degradation During Cyclic Loading[J].Granular Matter,2006, 8(3-4): 195-204.

[18] INDRARATNA B, THAKUR P, VINOD J. Experimental and Numerical Study of Railway Ballast Behavior under Cyclic Loading[J]. International Journal of Geomechanics, 2010, 10(4): 136-144.

[19] WILLIAMS J R, PENTL A P. Superquadrics and Modal Dynamics for Discrete Elements in Interactive Design [J]. Engineering Computations, 1992, 9(2): 115-127.

[20] CLEARY P W. Large Scale Industrial DEM Modeling[J]. Engineering Computations, 2004, 21(2/3/4): 169-204.

[21] HOHNER D, WIRTZ S, SCHERER V. A Numerical Study on the Influence of Particle Shape on Hopper Discharge within the Polyhedral and Muti-sphere Discrete Element Method[J]. Powder Technology, 2012, (226): 16-28.

[22] PENA A A, LIZCANO A, ALONSO-MARROQUIN F, et al. Biaxial Test Simulations Using a packing of Polygonal Particles[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(2):143-160.

[23] HOSSEININIA E S. Discrete Element Modeling of Inherently Anisotropic Granular Assemblies with Polygonal Particles[J]. Particuology, 2012,(10):542-552.

[24] GALINDO-TORRES S A, MUNOZ J D. Minkowski - Voronoi Diagrams as a Method to Generate Random Packings of Spheropolygons for the Simulation of Soils[J]. Physical Review E, 2010, 82(5): 1-12.

[25] ALONSO-MARROQUIN F, WANG Y C. An Efficient Algorithm for Granular Dynamics Simulations with Complex-shaped Objects [J]. Granular Matter, 2009, 11(5), 317-329.

[26] GALINDO-TORRES S A, ALONSO-MARROQUIN F, WANG Y C, et al. Molecular Dynamics Simulation of Complex Particles in Three Dimensions and the Study of Friction Due to Nonconvexity[J]. Physical Review E, 2009, 79(6): 060301.

[27] HOPKINS M A. Discrete Element Modeling with Dilated Particles[J]. Engineering Computations, 2003, 21(2/3/4): 422-430.

[28] POURNIN L, WEBER M, TSUKAHARA M,et al. Three-dimensional Distinct Element Simulation of Spherocylinder Crystallization[J]. Granular Matter, 2005, 7(2): 119-126.

[29] GALINDO-TORRES S A, PEDROSO D M, WILLIAMS D J,et al. Breaking Processes in Three-dimensional Bonded Granular Materials with General Shapes [J]. Computer Physics Communications, 2012, 183(2):266-277.

[30] YAN Y, JI S. Discrete Element Modeling of Direct Shear Tests for a Granular Material[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2010,34(9):978-990.

[31] JI S, SHEN H H. Effect of Contact Force Models on Granular Flow Dynamics[J]. Journal of Engineering Mechanics, 2006, 132(11): 1252-1259.

[32] JOHNSON K L. Contact Mechanics[M].Cambridge: Cambridge University Press，1987.

[33] 孫其誠, 王光謙. 顆粒物質(zhì)力學(xué)導(dǎo)論[M]. 北京: 科學(xué)出版社. 2009.