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

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

(1.廣西科技大學土木建筑工程學院，廣西柳州545006;2.同濟大學地下建筑與工程系，上海200092;3.廣西大學 工程防災與結構安全教育部重點實驗室，廣西南寧530004)

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

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

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

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

1.1 土工格柵的顆粒流模型

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

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

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

1.2 土樣細觀參數確定

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

圖2 數值試樣與標準砂顆粒級配對比Fig.2 Comparison of grade curves between numerical sample and standard sand

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

通過不斷調整kn、ks及fc的等細觀參數使數值試樣的宏觀力學表現逐漸逼近實際標準砂結果，具體數值為:顆粒密度2 650 kg/m3，最大粒徑3.75 mm，最小粒徑1.0 mm，孔隙率0.35，法向剛度2 ×106N/m，切向剛度2×106N/m，摩擦系數0.5。計算所得內摩擦角為29°，較接近標準砂室內試驗的內摩擦角27.3°。

圖3 數值三軸試驗試樣Fig.3 Triaxial test of PFC3D

1.3 筋土界面模型數值模擬步驟

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

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

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

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

圖4 土工合成材料拉拔試驗數值模型示意圖Fig.4 PFC3D model of geosynthetics pullout test

表1 模型箱各分區顆粒參數Table 1 Particle parameters of each partition model

2 宏觀響應分析

2.1 與室內模型試驗對比驗證

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

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

2.2 界面區域應力場分析

圖6、7分別給出了數值試驗拉拔模型中部XZ剖面土體中X向和Z向應力等值線云圖。

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

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

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

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

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

3 筋土界面細觀響應分析

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

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

圖8 顏色與位移對應關系圖Fig.8 The relationship map of color and displacement

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

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

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

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

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

圖11 室內拉拔試驗形成穩定界面的X向位移場Fig.11 The forming stability X displacement field of reinforced soil interface in laboratory pullout test

3.2 筋土界面顆粒細觀參數變化規律

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

所要量測的細觀參數為:局部孔隙率、接觸數和滑動比。其中局部孔隙率是指測量球中土體孔隙與整個球體體積的比值。局部孔隙率的變化是土體發生剪脹與剪縮最直接的反應。接觸數是測量球內單個顆粒與周圍顆粒接觸點數的平均值。實際上，顆粒間的接觸數主要是由于拉拔過程中被格柵帶動的顆粒與周圍顆粒相互擠壓咬合產生的，所以接觸數的變化規律一定程度上可反映拉拔過程中顆粒的運動特性。滑動比的定義為在測量球內發生滑動的接觸點數與總點數的比值。研究滑動比的變化規律可以進一步判別拉拔過程中顆粒與界面接觸發生位移和旋轉的劇烈程度。

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

3.2.1 靠近筋土界面區域細觀參數變化規律

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

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

此外，靠近界面區域土體的滑動比，在拉拔過程中呈現劇烈波動，表明顆粒在格柵帶動下發生了明顯的滑動和翻轉。

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

3.2.2 遠離界面區域土體的細觀參數變化規律

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

圖14 拉拔過程中遠離筋土界面區域細觀參數變化Fig.14 The change of mesoscopic parameters away from the reinforced soil interface

此外，對于接觸數變化情況，可看出遠離界面區域土體接觸數要遠小于界面區域土體接觸數，隨著拉拔位移增加，模型各區域的平均接觸數逐漸收斂于穩定的數值，表明顆粒間相互接觸最后達到較為均衡的狀態，整個系統趨于穩定狀態。

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

4 筋土界面荷載傳遞機理討論

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

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

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

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

5 結論

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

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

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

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

［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］夏才初，宋英龍，唐志成，等.粗糙節理剪切性質的顆粒流數值模擬［J］.巖石力學與工程學報，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］周健，史旦達，賈敏才，等.砂土單調剪切力學性狀的顆粒流模擬［J］.同濟大學學報:自然科學版，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.