加筋土拉拔界面作用的離散元細(xì)觀模擬

王家全，周健，吳輝琴，徐華

(1.廣西科技大學(xué)土木建筑工程學(xué)院，廣西柳州545006;2.同濟(jì)大學(xué)地下建筑與工程系，上海200092;3.廣西大學(xué) 工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，廣西南寧530004)

加筋土的界面摩擦特性直接影響加筋工程的安全與穩(wěn)定性，對其接觸界面進(jìn)行研究具有重要意義，有利于弄清土與土工合成材料之間相互作用的機(jī)理。數(shù)值仿真作為對室內(nèi)試驗(yàn)的有力補(bǔ)充，目前已被廣泛用于包括拉拔試驗(yàn)在內(nèi)的界面摩擦的仿真分析。國外較多學(xué)者［1-3］利用有限元方法，分別采用非線性硬化本構(gòu)模型和彈塑性模型研究筋土界面的相互作用，分析筋材的拉力、應(yīng)變等特征。有限元法可以較好的分析筋土界面的應(yīng)力應(yīng)變關(guān)系，但無法從顆粒層次上揭示界面的抗剪強(qiáng)度和筋土相互作用的細(xì)觀力學(xué)機(jī)理。

相比較于有限元方法等連續(xù)介質(zhì)力學(xué)的數(shù)值方法，應(yīng)用離散元方法對加筋土進(jìn)行分析的研究才剛剛開始，這方面的成果相對較少［4-6］。還有部分學(xué)者采用了離散元法研究接觸界面的細(xì)觀性狀［7-10］，主要是針對土與結(jié)構(gòu)界面的研究，從二維角度研究細(xì)觀參數(shù)變化對宏觀力學(xué)反應(yīng)的影響。

本文對顆粒離散元PFC3D進(jìn)行二次開發(fā)，建立模擬土工格柵拉拔試驗(yàn)的三維顆粒流細(xì)觀模型，研究筋土界面作用內(nèi)在機(jī)理，分析界面細(xì)觀參數(shù)的演化規(guī)律，從而更深入地了解拉拔過程中筋土界面宏觀力學(xué)響應(yīng)與細(xì)觀參數(shù)演化，有助于進(jìn)一步揭示筋土界面相互作用的內(nèi)在機(jī)理。

1 筋土界面作用三維顆粒流細(xì)觀模擬

1.1 土工格柵的顆粒流模型

在土工格柵拉拔試驗(yàn)的數(shù)值模型中，土工格柵采用2種不同粒徑的顆粒通過平行粘結(jié)模型組合而成，其中格柵的節(jié)點(diǎn)由粒徑較大的顆粒模擬，橫肋與縱肋采用粒徑較小的顆粒。圖1為拉拔界面試驗(yàn)采用的土工格柵離散元數(shù)值模型。

圖1 土工格柵離散元數(shù)值模型Fig.1 Geogrid of PFC3D

土工格柵在PFC3D模型中采用平行粘結(jié)模擬承受拉力，具體參數(shù)為:顆粒密度1 000 kg/m3，格柵粒徑5.0 mm，節(jié)點(diǎn)粒徑2.5 mm，法向剛度2 ×106N/m，切向剛度2×106N/m，法向粘結(jié)剛度9.9×1012N/m3，切向粘結(jié)剛度9.9 ×1012N/m3，摩擦系數(shù)0.5。所采用的細(xì)觀參數(shù)能基本反映土工格柵的受拉性狀。

1.2 土樣細(xì)觀參數(shù)確定

在數(shù)值模擬中，通常采用將實(shí)際砂土平均粒徑適當(dāng)放大的辦法來減少顆粒的數(shù)量，同時(shí)通過對最大、最小粒徑比的控制，使數(shù)值試樣的級配(不均勻系數(shù)Cu，曲率系數(shù)Cc)接近于實(shí)際標(biāo)準(zhǔn)砂［11］。室內(nèi)試驗(yàn)所用的標(biāo)準(zhǔn)砂是一種級配均勻的中砂，其平均粒徑d50=0.35 mm，不均勻系數(shù)Cu=1.560，曲率系數(shù)Cc=0.986。圖2為數(shù)值試樣與標(biāo)準(zhǔn)砂的級配曲線對比圖。由圖2可知，數(shù)值試樣d50=2.1 mm，Cu=1.533，Cc=0.936;與標(biāo)準(zhǔn)砂的級配相比，除d50增大外，Cu、Cc值均與標(biāo)準(zhǔn)砂相接近。

圖2 數(shù)值試樣與標(biāo)準(zhǔn)砂顆粒級配對比Fig.2 Comparison of grade curves between numerical sample and standard sand

本文采取的數(shù)值三軸試驗(yàn)來獲取與室內(nèi)試驗(yàn)性質(zhì)接近的砂土顆粒細(xì)觀參數(shù)，如圖3所示;利用伺服機(jī)制施加所需圍壓，然后施加軸向位移直到試樣破壞。

通過不斷調(diào)整kn、ks及fc的等細(xì)觀參數(shù)使數(shù)值試樣的宏觀力學(xué)表現(xiàn)逐漸逼近實(shí)際標(biāo)準(zhǔn)砂結(jié)果，具體數(shù)值為:顆粒密度2 650 kg/m3，最大粒徑3.75 mm，最小粒徑1.0 mm，孔隙率0.35，法向剛度2 ×106N/m，切向剛度2×106N/m，摩擦系數(shù)0.5。計(jì)算所得內(nèi)摩擦角為29°，較接近標(biāo)準(zhǔn)砂室內(nèi)試驗(yàn)的內(nèi)摩擦角27.3°。

圖3 數(shù)值三軸試驗(yàn)試樣Fig.3 Triaxial test of PFC3D

1.3 筋土界面模型數(shù)值模擬步驟

拉拔數(shù)值模型采用與室內(nèi)模型相同的尺寸，即0.6 m ×0.4 m ×0.4 m(X×Y×Z)。

1)砂土顆粒隨機(jī)生成于由六面墻圍成的長方形區(qū)域;考慮當(dāng)前計(jì)算機(jī)水平，采用分塊建模原則，如圖4所示，在中心界面區(qū)域1中生成粒徑較小的顆粒，面積為X×Z=0.6 m ×0.1 m，在區(qū)域2中生成粒徑較大的顆粒，面積為X×Z=0.6 m×0.1 m，在區(qū)域3生成粒徑較大的顆粒，面積為X×Z=0.6 m×0.1 m。模型各區(qū)顆粒參數(shù)詳見表1，砂土的細(xì)觀參數(shù)取自1.2節(jié)的三軸試驗(yàn)擬定的參數(shù)。

2)在模型中間位置生成土工格柵數(shù)值模型，土工格柵沿X方向水平設(shè)置，土工合成材料埋入模型箱內(nèi)0.55 m，并利用伺服機(jī)制施加指定的壓力。

3)將模型所有顆粒的位移和速率清零，對土工格柵施加X向恒定速率1 mm/min，將土工格柵拔出25 mm，并記錄各種變量。

圖4 土工合成材料拉拔試驗(yàn)數(shù)值模型示意圖Fig.4 PFC3D model of geosynthetics pullout test

表1 模型箱各分區(qū)顆粒參數(shù)Table 1 Particle parameters of each partition model

2 宏觀響應(yīng)分析

2.1 與室內(nèi)模型試驗(yàn)對比驗(yàn)證

為驗(yàn)證數(shù)模模擬的可靠性，與室內(nèi)模型做對比，2種模型采用相同尺寸和法向壓力均為30 kPa，土工格柵拉拔速度均為1 mm/min，圖5中給出了數(shù)值試驗(yàn)與室內(nèi)試驗(yàn)在法向壓力30 kPa下拉拔阻力與位移的對比曲線。由圖5可知，室內(nèi)試驗(yàn)土工格柵拉拔力隨位移增長而增長，到達(dá)峰值后，開始緩慢回落，數(shù)值試驗(yàn)也呈現(xiàn)出相同的規(guī)律，且二者拉拔阻力峰值較為接近，顆粒流(particle flow code，PFC)數(shù)值細(xì)觀模型試驗(yàn)較好的重現(xiàn)了室內(nèi)試驗(yàn)受力結(jié)果。巖土體宏觀力學(xué)特征規(guī)律是細(xì)觀參數(shù)綜合響應(yīng)的結(jié)果，說明所建立的PFC細(xì)觀數(shù)值模型能有效模擬室內(nèi)拉拔試驗(yàn)過程。

圖5 土工格柵的拉拔阻力－位移曲線Fig.5 The curves of geogrid pullout load-displacement

2.2 界面區(qū)域應(yīng)力場分析

圖6、7分別給出了數(shù)值試驗(yàn)拉拔模型中部XZ剖面土體中X向和Z向應(yīng)力等值線云圖。

通過圖6、7可以看出，在土工格柵拉拔開始前，X向、Z向應(yīng)力在模型中部界面上下區(qū)域分布完全不同，界面上部應(yīng)力等值線較為稀疏，而下部等值線密集，表明由于自重的原因，下部顆粒相對上部而言承受較大應(yīng)力，符合實(shí)際情況。

隨著拉拔的進(jìn)行，在筋土界面上下4 cm范圍內(nèi)形成了X向應(yīng)力集中區(qū)域，該應(yīng)力集中區(qū)為土工格柵表面摩擦力的主要來源，此后，隨著拉拔位移S進(jìn)一步增大，該X向應(yīng)力集中區(qū)往前平移，靠近前墻的應(yīng)力進(jìn)一步增大，而靠近后墻應(yīng)力逐漸較小，主要是由于土工格柵拔出后，造成格柵后部脫空，顆粒跌落補(bǔ)充空隙，造成相應(yīng)區(qū)域應(yīng)力減小。

此外，從圖7可以看出，在拉拔過程中Z向應(yīng)力的變化規(guī)律略有不同，隨著拉拔位移增大，在沿土工格柵界面區(qū)域出現(xiàn)了幾處Z向應(yīng)力集中區(qū)，說明土工格柵橫肋的存在，提供了較顯著的被動承載力，拉拔過程橫肋擠壓土體，造成局部區(qū)域應(yīng)力增大，橫肋的被動承載力逐漸增大。

圖6 拉拔過程模型箱XZ剖面處X向應(yīng)力等值線Fig.6 The X stress field of model box XZ section

圖7 拉拔過程模型箱XZ剖面處Z向應(yīng)力等值線Fig.7 The Z stress field of model box XZ section

3 筋土界面細(xì)觀響應(yīng)分析

3.1 筋土界面顆粒位移演化分析

在離散元PFC中，采用FISH語言開發(fā)位移云圖顯示程序，位移的大小通過不同顆粒顏色來進(jìn)行表示，獲知模型任意區(qū)域顆粒的位移分布情況。圖8為顏色與位移量大小的關(guān)系對應(yīng)圖。

圖8 顏色與位移對應(yīng)關(guān)系圖Fig.8 The relationship map of color and displacement

圖9為PFC數(shù)值模型箱中部XZ剖面砂土X向位移演化圖。在拉拔開始階段，筋土界面區(qū)域較大范圍內(nèi)的顆粒均發(fā)生了位移，當(dāng)拉拔位移達(dá)到10 mm時(shí)，界面區(qū)域顆粒位移基本上發(fā)展為一個(gè)狹長的位移集中帶，拉拔位移達(dá)到20 mm時(shí)，該位移集中帶已趨于穩(wěn)定，基本呈上下對稱形狀。

數(shù)值模擬與室內(nèi)試驗(yàn)相比，可以得到模型任何位置顆粒的位移情況。圖10為拉拔20 mm時(shí)模型箱中部YZ剖面的筋土界面顆粒位移場，在XZ剖面的筋土界面處亦形成了一個(gè)均勻的位移集中帶，其大小與圖9中XZ剖面觀察到的集中帶基本一致，表明土工格柵在拉拔過程中，帶動了整個(gè)筋土界面處砂土的位移，形成了穩(wěn)定的位移集中帶。

圖11為室內(nèi)拉拔試驗(yàn)形成穩(wěn)定界面的位移場，對比圖9、10離散元數(shù)值模擬結(jié)果，室內(nèi)試驗(yàn)和數(shù)值試驗(yàn)均形成了相同的位移集中帶，從另一側(cè)面驗(yàn)證數(shù)值模型的正確。數(shù)值試驗(yàn)中格柵上界面位移帶厚度略大于下界面位移帶厚度，這與室內(nèi)試驗(yàn)結(jié)果較為一致，主要原因是顆粒在自重和法向壓力作用下，下部界面顆粒較為密實(shí)，拉拔所能帶動的顆粒范圍要略小于上部界面區(qū)域。

圖9 模型箱XZ剖面處X向位移場演化圖Fig.9 The X displacement field of model box XZ section

圖10 模型箱YZ剖面處形成穩(wěn)定筋土界面的位移場Fig.10 The forming stability X displacement field of reinforced soil interface in model box YZ section

圖11 室內(nèi)拉拔試驗(yàn)形成穩(wěn)定界面的X向位移場Fig.11 The forming stability X displacement field of reinforced soil interface in laboratory pullout test

3.2 筋土界面顆粒細(xì)觀參數(shù)變化規(guī)律

在室內(nèi)試驗(yàn)中由于量測手段的限制，一直難以精確描述模型內(nèi)部筋土界面作用規(guī)律，顆粒流三維數(shù)值模型可以精確的獲取拉拔過程中任意區(qū)域的細(xì)觀參數(shù)變化。為量測拉拔過程中界面土體細(xì)觀參數(shù)演化規(guī)律，在模型箱中部XZ剖面處設(shè)置8排測量球，圖12為測量球布置示意圖，其中靠近筋土界面區(qū)域測量球編號為3、4、5、6，遠(yuǎn)離界面區(qū)域測量球?yàn)?、2、7、8，測量球直徑為2 cm。將每一排2個(gè)測量球數(shù)據(jù)平均作為該排的細(xì)觀參數(shù)。

所要量測的細(xì)觀參數(shù)為:局部孔隙率、接觸數(shù)和滑動比。其中局部孔隙率是指測量球中土體孔隙與整個(gè)球體體積的比值。局部孔隙率的變化是土體發(fā)生剪脹與剪縮最直接的反應(yīng)。接觸數(shù)是測量球內(nèi)單個(gè)顆粒與周圍顆粒接觸點(diǎn)數(shù)的平均值。實(shí)際上，顆粒間的接觸數(shù)主要是由于拉拔過程中被格柵帶動的顆粒與周圍顆粒相互擠壓咬合產(chǎn)生的，所以接觸數(shù)的變化規(guī)律一定程度上可反映拉拔過程中顆粒的運(yùn)動特性?；瑒颖鹊亩x為在測量球內(nèi)發(fā)生滑動的接觸點(diǎn)數(shù)與總點(diǎn)數(shù)的比值。研究滑動比的變化規(guī)律可以進(jìn)一步判別拉拔過程中顆粒與界面接觸發(fā)生位移和旋轉(zhuǎn)的劇烈程度。

圖12 量測細(xì)觀參數(shù)的測量球布置示意圖Fig.12 Measurement ball layout map for measuring mesoscopic parameters

3.2.1 靠近筋土界面區(qū)域細(xì)觀參數(shù)變化規(guī)律

圖13為靠近筋土界面區(qū)域土體的細(xì)觀參數(shù)變化曲線。從圖13可以發(fā)現(xiàn)，測量球3和測量球4分別位于筋土上下界面處，其孔隙率隨拉拔位移增長而增大，其中在測量球3區(qū)域拉拔位移約6.5 mm時(shí)，孔隙率增大到峰值，測量球4區(qū)域在拉拔位移20 mm時(shí)孔隙率達(dá)到峰值，然后開始回落。由此可以發(fā)現(xiàn)，在上下界面區(qū)域均存在局部孔隙率隨位移增大的情況，界面局部位置發(fā)生剪脹，結(jié)合3.1節(jié)的拉拔阻力與拉拔位移變化曲線(圖5)，拉拔阻力峰值約處于拉拔位移5.2 mm位置，兩者峰值對應(yīng)位移較為接近，表明孔隙率增大，界面局部區(qū)域土體發(fā)生剪脹應(yīng)力增大，促使拉拔阻力增大，而后孔隙率慢慢減小，拉拔阻力回落。

對應(yīng)接觸數(shù)變化情況，測量球3和測量球4的接觸數(shù)與局部孔隙率剛好相反，隨拉拔位移增大而減小，接觸數(shù)先是減小到谷值，然后慢慢回升。區(qū)域顆粒隨拉拔發(fā)生剪脹而后，從而造成接觸數(shù)減少，此后孔隙率減小，接觸數(shù)相應(yīng)得到恢復(fù)。模擬結(jié)果還可以看出，模型下部(如測量球1、2區(qū)域)接觸數(shù)變化甚微，而上部(如測量球7、8區(qū)域)接觸數(shù)開始階段變化略大，但均遠(yuǎn)小于界面區(qū)域土體接觸數(shù)變化情況。

此外，靠近界面區(qū)域土體的滑動比，在拉拔過程中呈現(xiàn)劇烈波動，表明顆粒在格柵帶動下發(fā)生了明顯的滑動和翻轉(zhuǎn)。

圖13 拉拔過程中靠近筋土界面區(qū)域細(xì)觀參數(shù)變化Fig.13 The change of mesoscopic parameters near the reinforced soil interface

3.2.2 遠(yuǎn)離界面區(qū)域土體的細(xì)觀參數(shù)變化規(guī)律

圖14為遠(yuǎn)離筋土界面區(qū)域細(xì)觀參數(shù)變化曲線，從圖中可以發(fā)現(xiàn)，遠(yuǎn)離筋土下界面的測量球1、2區(qū)域，孔隙率波動很小，表明格柵拉拔對該區(qū)域土體影響甚微。而遠(yuǎn)離上界面的測量球7、8區(qū)域，其孔隙率變化要比下界面對應(yīng)位移的測量球1、2區(qū)域的孔隙率要稍大，先是減少，而后開始緩慢上升保持平穩(wěn)趨勢。由此可以發(fā)現(xiàn)格柵拉拔造成的上界面影響區(qū)域要大于下界面的影響區(qū)域。

圖14 拉拔過程中遠(yuǎn)離筋土界面區(qū)域細(xì)觀參數(shù)變化Fig.14 The change of mesoscopic parameters away from the reinforced soil interface

此外，對于接觸數(shù)變化情況，可看出遠(yuǎn)離界面區(qū)域土體接觸數(shù)要遠(yuǎn)小于界面區(qū)域土體接觸數(shù)，隨著拉拔位移增加，模型各區(qū)域的平均接觸數(shù)逐漸收斂于穩(wěn)定的數(shù)值，表明顆粒間相互接觸最后達(dá)到較為均衡的狀態(tài)，整個(gè)系統(tǒng)趨于穩(wěn)定狀態(tài)。

與界面區(qū)域情況相比，遠(yuǎn)離界面區(qū)域的土體滑動比波動很小，表明越遠(yuǎn)離界面區(qū)域，滑動比變化越小，其中模型下部變化要大大小于上部變化，從另一側(cè)面反映出模型下部顆粒受上部土體重力壓密作用，幾乎不發(fā)生位移和旋轉(zhuǎn)。

4 筋土界面荷載傳遞機(jī)理討論

拉拔試驗(yàn)不僅給工程設(shè)計(jì)提供了界面摩擦系數(shù)和極限拉拔力，并且有助于深入研究筋土界面作用荷載傳遞的機(jī)理。

土工格柵拉拔阻力PR主要由摩擦阻力和橫肋提供的被動承載阻力組成:

式中:PRS為格柵表面摩阻力，PRS=(csg+σ'tan φsg)A，csg為界面粘聚力，σ＇為有效法向應(yīng)力，φsg為界面摩擦角，A為界面總面積;PRB為格柵橫肋的被動承載力。

結(jié)合前面幾節(jié)分析，隨著拉拔位移的逐漸增加，土工格柵表面的摩擦阻力以及橫肋被動承載力隨之逐步發(fā)揮出來，特別是筋土界面處X向應(yīng)力高度集中，當(dāng)拉拔位移達(dá)到一定程度，格柵表面摩擦阻力PRS基本保持不變，隨著拉拔位移的進(jìn)一步增大，格柵橫肋繼續(xù)擠壓土體，格柵橫肋被動承載力PRB開始發(fā)揮關(guān)鍵作用，在界面橫肋擠壓區(qū)域形成了若干個(gè)Z向應(yīng)力集中區(qū)，該階段拉拔阻力的增長主要依靠橫肋的被動承載力來提供，在細(xì)觀上，橫肋擠壓剪切土體，細(xì)觀參數(shù)分析中發(fā)現(xiàn)界面區(qū)孔隙率增大，土體內(nèi)部剪脹，拉拔阻力峰值與界面區(qū)孔隙率峰值較為接近，在宏觀力學(xué)數(shù)值與細(xì)觀孔隙數(shù)值的演化上存在著對應(yīng)關(guān)系;當(dāng)格柵的表面摩擦阻力和橫肋被動承載力達(dá)到極限狀態(tài)時(shí)，筋土界面周圍土體發(fā)生剪切破壞，格柵以均勻的速度被拔出。

5 結(jié)論

1)離散元PFC3D經(jīng)過二次開發(fā)，可有效模擬土工合成材料與散體顆粒的拉拔試驗(yàn)，深入探討加筋土的細(xì)觀機(jī)理。

2)在筋土界面狹長區(qū)域，砂土主要以水平位移為主，在遠(yuǎn)離界面區(qū)域，砂土以法向位移為主。在土工格柵附近形成了上下兩個(gè)界面位移集中帶，上界面影響厚度略大于下界面影響厚度。

3)在筋土界面附近區(qū)域形成了高度的X向應(yīng)力集中區(qū)，而沿格柵長度方向處界面區(qū)域的Z向應(yīng)力出現(xiàn)幾個(gè)應(yīng)力集中區(qū)，表明對應(yīng)部位格柵橫肋提供的被動承載力較為顯著。

4)在筋土界面區(qū)域，細(xì)觀上局部孔隙率隨拉拔位移增大而增大，接觸數(shù)減少，滑動比上下波動劇烈，界面局部位置發(fā)生剪脹，宏觀上對應(yīng)拉拔阻力增大，而后界面區(qū)域孔隙率慢慢減小，拉拔阻力回落。

［1］YOGARAJAH I，YEO K C.Finite element modeling of pull-out tests with load and strain measurements［J］.Geotextiles and Geomembranes，1994，13(1):43-54.

［2］MOSTAFA A，EL S.Behavior of strip footing on geogrid-reinforced sand over a soft clay slope［J］.Geotextiles and Geomembranes，2007，25(1):50-60.

［3］BERGADO D T，TEERAWATTANASUK C.2D and 3D numerical simulations of reinforced embankments on soft ground［J］.Geotextiles and Geomembranes，2008，26(1):39-55.

［4］RICHARD P J，PETER J B，MICHAEL E P，et al.DEM Simulation of granular media-structure interface:effects of surface roughness and particle shape［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1999，23(6):531-547.

［5］CALVETTI F，PRISCO D C，NOVA R.Experimental and numerical analysis of soil-pipe interaction［J］.Journal of Geotechnical and Geoenvironmental Engineering，2004，130(12):1292-1299.

［6］夏才初，宋英龍，唐志成，等.粗糙節(jié)理剪切性質(zhì)的顆粒流數(shù)值模擬［J］.巖石力學(xué)與工程學(xué)報(bào)，2012，31(8):1545-1552.XIA Caichu，SONG Yinglong，TANG Zhicheng，et al.Particle flow numerical simulation for shear behavior of rough joints［J］.Chinese Journal of Rock Mechanics and Engineering，2012，31(8):1545-1552.

［7］HIROAKI T，MASATOSHI M，AKIRA O，et al.Simulation of soil deformation and resistance at bar penetration by the distinct element method［J］.Journal of Terramechanics，2000，37(1):41-56.

［8］JIANG M J，LEROUEIL S，KONRAD J M.Yielding of microstructured geomaterial by distinct element method analysis［J］.Journal of Engineering Mechanics，2005，131(11):1209-1213.

［9］CAMBORDE F，MARIOTTI C，DONZE F V.Numerical study of rock and concrete behaviour by discrete element modelling［J］.Computers and Geotechnics，2000，27(4):225-247.

［10］JENSEN R P，EDIL T B，BOSSCHER P J，et al.Effect of particle shape on interface behavior of DEM-simulated granular materials［J］.International Journal of Geomechanics，2001，1(1):1-19.

［11］周健，史旦達(dá)，賈敏才，等.砂土單調(diào)剪切力學(xué)性狀的顆粒流模擬［J］.同濟(jì)大學(xué)學(xué)報(bào):自然科學(xué)版，2007，35(10):1299-1304.ZHOU Jian，SHI Danda，JIA Mincai，et al.Numerical simulation of mechanical response on sand under monotonic loading by particle flow code［J］.Journal of Tongji University:Natural Science，2007，35(10):1299-1304.