加筋土拉拔界面作用的离散元细观模拟

王家全，周健，吴辉琴，徐华

(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.