基于土工格栅—砾石室内大型直剪试验及其结果，利用离散元软件，模拟法向循环荷载下的筋土界面直剪试验过程。试样尺寸与直剪仪剪切盒相同，上、下盒分别为800 mm×200 mm×50 mm和600 mm×200 mm×50 mm。本次建立颗粒流模型时，采用试验实际用土的级配进行土颗粒的模拟，颗粒参数（法向刚度、切向刚度和摩擦系数等）则采用“试错标定法”^[17]，通过一系列数值试验并将模拟结果与室内试验结果进行对比，经过不断调整、试算、拟合后，最终得到比较吻合宏观力学结果的参数取值。粗粒土颗粒采用线性接触模型，颗粒采用球颗粒，格栅单元采用平行黏结接触模型，如图1所示。将标定的细观参数分别赋予粗粒土颗粒和格栅颗粒，具体参数取值如表1所示。

图1 土工格栅示意图Fig.1 Schematic diagram of geogrid

利用离散元软件进行模拟时，墙体只能施加位移和速度，无法施加荷载，因此，法向循环荷载不能通过对墙施加的方式实现。为了实现法向荷载的施加，通过生成颗粒组成的clump刚性簇面板，对其中颗粒施加循环荷载，以达到对试样施加法向循环荷载的目的。借助离散元软件中的FISH编程功能，使试样受到的动荷载与室内试验相同。

表1 数值模拟参数Table 1 Parameters of numerical simulation

模拟直剪试验的过程主要有以下几个步骤：

（1）在剪切盒尺寸的空间内，按级配生成随机的球形颗粒，并标定合适的参数，在剪切面上和两侧分别生成测量球，如图2所示。

图2 测量球位置Fig.2 Position of measuring ball

（2）在上、下剪切盒界面处生成土工格栅模型，并标定合适的参数。

（3）在自重状态下进行初始状态的平衡运算。

（4）生成clump面板，待初始荷载稳定后，对其施加法向循环荷载，模型如图3所示。

图3 直剪试验数值模型Fig.3 Numerical model of direct shear tests

（6）通过控制下盒向右作水平运动，实现剪切过程。

为验证三维离散元数值模型的准确性，将室内循环荷载（40~80 kPa）下的筋土界面直剪试验结果与数值模拟结果进行对比，结果如图4所示。可见，剪切过程中剪切应力与体变特征均与试验高度重合，证明了数值模型的正确性。

为研究法向循环荷载作用下筋土界面的变形特征，需要关注界面局部化、非连续的剪切带演化。因此，本文在剪切过程开始前，设置12列竖向宽度相同的黑色颗粒带，剪切过程各列的颗粒演化如图5所示。

在剪切位移为0.1 mm时，即剪切初期时，竖向条带基本没有变形，处于纯剪切主导的状态。剪切应力到达峰值附近时（剪切位移为5 mm），竖向条带的变形也较小，只在侧边的竖列发生轻微的变化。在峰后软化阶段（剪切位移30 mm），局部变形现象愈发明显，从侧边至中间的竖向条带均演化出应变局部化现象，体现出此时筋土界面非线性变化的特征。剪切结束后（剪切位移60 mm），在剪切面上、下两侧形成了局部变形集中的剪切带，竖向颗粒带的弯曲明显，即图中虚线的区间范围。在法向循环荷载作用下的筋土界面直剪试验过程中，峰值前以纯剪切为主，变形较小，而峰后残余阶段以非线性应变为主导，局部变形较明显。

图4 数值模拟与室内试验Fig.4 Comparison of numerical simulation and indoor experimental results

图5 剪切过程剪切带演化Fig.5 Shear band evolution in shear proces

在颗粒物质力学中，颗粒体系内部的接触互相连接逐渐形成能够传递荷载的线性链状结构，即为力链。这些力链错综复杂地在体系内部形成力链网络，力链网络是土骨架的粒间接触力的宏观表现形式，可用于反映颗粒体系受外荷载作用时受力响应机制的变化。本文通过分析剪切过程中颗粒力链网络的演化，研究了土骨架的受力情况。

图6为法向循环荷载下界面剪切过程力链演化图。在完成初始应力加载、剪切初期时（剪切位移0.1 mm），由于竖向荷载和下盒剪切力的作用，下盒的颗粒力链部分强于上盒。经过一段时间剪切后，力链的强弱状态发生了改变，力链网络发生重构。当达到峰值阶段（剪切位移5 mm），上盒右侧力链相较于其他位置更强，而下盒则相反。当试验进入峰后软化阶段，上述力链状态继续强化，试样强力链进一步集中，主要在上盒右侧方、下盒左下和中间位置。从试样立体状态看（图7，剪切位移60 mm），在界面形成了贯通上、下剪切盒倾斜的强力链，这表明剪切位移大小和剪切方向会影响强力链的产生和倾斜角度。

图6 剪切过程力链演化Fig.6 Force chain evolution in shear process

图7 剪切完成时立体角度的力链状态Fig.7 Force chain state of solid angle with cutting completed

图8为界面剪切过程位移场变化图。在初始剪切阶段（剪切位移0.1 mm），内部位移场处于一种无序的状态，且存在大量向下位移的颗粒，对应着试样初始剪缩状态。随着剪切的继续，下盒向下运动的颗粒逐渐减少，而转向沿剪切速度方向向右运动。上盒的颗粒由起始阶段各个方向的运动状态均存在，渐渐旋转角度变为以向上和向右运动为主（剪切位移5 mm），这主要是下盒颗粒挤压和摩擦导致的，证明了持续剪胀现象的发生，而到了残余应力状态，又出现了部分向下运动的颗粒，对应着剪胀位移放缓的现象。

图8 剪切过程位移场演化Fig.8 Evolution of displacement field in shear process

由于剪切作用带动筋土界面颗粒产生了位移等原因，孔隙率一直处于变化中，量测区域的局部孔隙率是试样紧实程度的细观表征，也是筋土界面剪胀性的直接证明。图9为不同区域的孔隙率改变量随剪切位移变化的曲线。

在初始剪切阶段，试样孔隙率减小为主，代表试样初始短暂的剪缩状态；第二阶段主要以孔隙率增加为主，表明此时剪胀位移更明显；最后剪胀变化趋于稳定的状态，相应的部分位置出现孔隙率减小的情况。同时，从不同的测量圆位置看，剪切面（位置4、5、6）的孔隙率增长现象最显著，其次是上剪切盒（位置1、2、3），而下盒（位置7、8、9）的试样孔隙变化较小，甚至局部出现较长的减缩过程。法向循环荷载与水平剪切作用下，格栅界面处受影响程度最大，剪胀性最显著，孔隙率增大从剪切带部分逐渐向上、下两侧发展，试样剪胀性主要受剪切面局部区域的影响。

图9 剪切过程孔隙率—位移曲线Fig.9 Porosity-displacement curve in shear process

配位数指颗粒的平均接触数目，是颗粒流模型分析中常用的参数。配位数越大，颗粒间的受力传递越充分，其大小一定程度上反映了试样所处的状态（如应力水平）。土力学中剪胀现象的本质是颗粒体系在外荷载作用下发生重排列，而颗粒的滑移、转动必然会导致试样内配位数发生变化。

图 10为不同位置测量圆内监测的平均接触数。剪切过程中，配位数呈现动态的变化，且不同区域内的的配位数变化差异显著。以位置1为例，配位数随剪切位移增加先波动着降低，而后在配位数1.75附近稳定变化，但在其他位置配位数的变化规律则不完全相同。位置3是上剪切盒中配位数最高的，而这也是前述力链场中力链较强的部分。在剪切面位置（图 10（b）），左、右两侧（位置4、6）的配位数波动程度更小，约在4左右变化，而中间位置（位置5）的配位数随位移的变化上、下起伏较大，且整体的值也更大。在下剪切盒（图 10（c）），左侧配位数（位置7）主要表现出减小的趋势，中间（位置8）的配位数随位移的增加呈先增大后减小的特征，而位置9的配位数变化较小；左侧配位数最高，而右侧配位数最小，这也能很好对应力链场的特征。总的来说，上剪切盒配位数曲线的振荡程度较大，而下盒配位数振荡较小；上盒试样配位数整体上小于下盒。

图 10 剪切过程配位数—位移曲线Fig.10 Coordination number-displacement curve in shear process

图 11为位置2和8的法向静、动荷载下的孔隙率随位移变化的曲线。动荷载作用下的孔隙率变化处于相应上、下限的静荷载作用下孔隙率之间，且在同样的位置，静、动力下的孔隙率随位移增减变化的趋势接近，可见动荷载下的试样密实度状态处于上、下界静荷载间。

图 12为位置2和位置8的法向静、动荷载下的配位数变化的比较曲线。与孔隙率的情况不同，动荷载下的配位数曲线并没有完全在相应的上、下界静荷载下的配位数曲线之间，而是更偏向下界荷载的状态，甚至低于下界荷载的配位数。在同一位置，曲线的增减变化规律也不完全一致，尤其是上盒颗粒在静力和动力下的差异性更显著，而下盒配位数变化更稳定，变化趋势接近。这是由于法向荷载的频繁变化，使得颗粒运动相较于静荷载更剧烈，进而导致颗粒接触更加不稳定和复杂，导致动荷载下的界面剪切时颗粒接触数偏小。

图 11 静、动力荷载下孔隙率比较Fig.11 Comparison of porosity under static and dynamic loads

图 12 静、动力荷载下配位数比较Fig.12 Comparison of coordination numbers under static and dynamic loads