基于颗粒流的膨胀土裂隙填充细观试验研究

王少鹏，张文慧，张福海，白岩辉，陈 宇

（1.河海大学岩土力学与堤坝工程教育部重点实验室，南京 210098；2.河海大学江苏省岩土工程技术工程研究中心，南京 210098）

膨胀土是由蒙脱石和伊利石等亲水性黏土矿物组成的高塑性黏土，其与普通黏土的力学特性不同，常具有明显的胀缩性、多裂隙性、强度衰减性和风化性等不良特性［1-4］.燕守勋等［5］利用XRD试验研究了膨胀土中蒙脱石含量与膨胀土膨胀势指标的关系；戴张俊等［6］通过SEM试验对南水北调中线工程中膨胀土的微观结构进行了定性分析；王冲、廖雄华等［7-8］利用颗粒流软件对黏土的细观特性进行探究；张福海、杨升、荆鹏等［9-13］将室内试验与模拟相结合，确保了模拟的准确性；冯晓腊、尹小涛等［14-18］研究了颗粒粒径和尺寸效应的影响.大量的现场资料表明，膨胀土边坡中含有多种类型的裂隙，裂隙中灰白色物质的特性对膨胀土边坡的稳定性影响明显［19-20］.但在细观条件下，膨胀土内部的充填裂隙对膨胀土颗粒细观特征等规律尚无系统的研究.

本文基于颗粒流法，在室内试验的基础上，通过Pfc二维颗粒流程序建立蒙脱石充填裂隙膨胀土直接剪切模型.将室内直剪试验的结果与数值模拟结果进行系统的比较，在细观条件下分析裂隙填充厚度对颗粒的运动、力链等细观的变化规律.

1 室内试验样本

本次试验所用的膨胀土样本取自江苏省南京市某膨胀土边坡支护工程，对土样经过处理后进行基础物理力学性质试验，具体试验参数如表1所示.由于非饱和土样的强度指标远大于饱和土样，因此本次试验所用土样均处于饱和状态.土样中的充填裂隙材料采用钠基膨润土进行模拟，其化学成分如表2所示.试验用膨润土中的蒙脱石含量为85%.由现场观测资料可知，原状土的天然含水率为26%，充填裂隙的天然含水率为30%，因此按照30%来配置充填裂隙的含水率.

表1 膨胀土样基础物理性质指标Tab.1 Physical property index of expansive soil foundation

表2 膨润土化学成分Tab.2 Chemical composition of expansive soil 单位：%

本次试验分别以充填裂隙厚度和蒙脱石含量为自变量进行探究.实际工程中充填裂隙的厚度大多在0～2 mm，裂隙中的蒙脱石含量明显高于裂隙两侧的土体.故本次模拟的填充层厚度定为1.0、1.5、2 mm，蒙脱石含量定为50%.出于工程安全考虑，充填裂隙应置于最不利的位置，故将充填裂隙的位置放在剪切面附近.本次按照土工试验规程进行直剪试验，每个试验均在100 kPa围压下进行.

2 Pfc直剪试验计算模型

2.1 建立模型

本文在建立模型时采用Pfc（version 5.0）二维建模软件.实际直剪试验中环刀的尺寸61.8 mm×20 mm，为更好地反映实际剪切过程，将模拟时剪切盒的几何尺寸与室内试验剪切盒的尺寸保持一致.

膨胀土颗粒中黏粒含量较高，颗粒的实际粒径较小，因此在建模时如果采用实际的颗粒半径，将会生成百万量级的颗粒.故在考虑运算速度、计算精度等因素的情况下，确定墙体内部的膨胀土球形颗粒个数为14 203个，且试样中颗粒粒径满足ASTM D3080-90中对高、宽与最大颗粒粒径的比分别不小于6和10的要求［21］，生成的颗粒如图1 a所示.在建立模型时将剪切盒内部分为上、中、下三部分，上部和下部都为膨胀土颗粒，中部为蒙脱石充填裂隙的土样，蒙脱石填充层的形态如图1 b所示（以2 mm为例）.颗粒的运动通过墙体进行约束，墙体与颗粒间通过接触发生作用，对于每一个颗粒，都能满足运动方程，颗粒的运动示意图如图1 c所示.由于墙体没有质量，因此墙体受到的不平衡力无法影响墙体的运动.整个模型由8面墙体进行组成，为防止在颗粒剪切过程中冲出墙体，将墙体的参数取较大值，法向刚度kn为3×106，刚度ks为1.5×106，墙体的分布如图1 d所示.

图1 模型示意图Fig.1 Model diagram

2.2 细观参数标定及模型选取

试验用土为弱膨胀土，土中黏粒含量较高，因此可选用接触黏结模型.接触黏结模型属于点接触模型，不需在计算中给出黏结半径等参数，且没有考虑颗粒间接触的法向应力，模型具有存在和消失两种状态.一旦在剪切向或者法向超过强度极限，则颗粒间约束消失，转化为无约束状态，此时颗粒相互作用遵从滑移模型.颗粒流模型的参数无法通过室内土工试验直接获得，需要通过伺服控制和数值试验对细观参数不断调整，膨胀土具体的细观模拟参数见表3.中间部分的蒙脱石填充区域颗粒参数与膨胀土颗粒区域参数不同，单独设置.

表3 膨胀土颗粒细观参数Tab.3 Mesoscopic parameters of expansive soil particles

孔隙比作为土体的重要参数，应当在建模时予以考虑.本文膨胀土的孔隙比参数的选取是依据何咏睿［22-23］等研究成果，将三维条件下的孔隙率转换为二维条件下的孔隙率，转换后的孔隙率为0.195.

3 试验结果分析

3.1 颗粒流试验结果标定

本次试验均是在竖向荷载为100 kPa下进行，以室内试验中的蒙脱石裂隙填充厚度的应力应变曲线为标定依据，对数值试验的结果进行调整，利用FISH函数记录每个时间步长的位移剪应力，并将实际值与模拟值进行比较，比较结果如图2所示.

当填充层厚度为1 mm时，室内试验得到的应力应变曲线在5 mm附近出现峰值，且残余强度值明显，而数值模型模拟的应力应变曲线在应变达到5 mm应力值在42.5 kPa附近上下波动，且残余强度不明显；当填充层厚度为1.5 mm时，模拟试验的峰值强度明显大于室内试验的峰值，模拟试验的残余强度与同等位移下的模拟值基本相等；当填充层厚度为2 mm时，模拟试验的峰值略大于室内试验的峰值，模拟试验的残余强度与同等位移下的模拟值相似.数值模拟的剪应力最大值明显小于室内试验的最大值，但整体上来看，数值模型与室内试验应力随应变的变化趋势大致相同，且最大值相差较小，因此模拟试验的效果理想.

图2 室内试验与模拟试验应力-应变对比图Fig.2 Stress-strain comparison between laboratory test and simulation test

3.2 颗粒位移分析

由图2应力应变对比可以看出，颗粒在剪切过程中根据坚硬应力的变化可分为三个阶段：第一阶段为剪应力增加阶段，第二阶段为剪应力峰值阶段，第三阶段为剪应力减小阶段，因此剪切位移取2、5、15 mm进行颗粒三个阶段的位移分析.图3展示了不同填充厚度下，随着剪切位移的变化土体内部单个颗粒位移的等值线云图变化情况.由图可知，在剪切第一阶段，不同的填充层厚度对颗粒位移的影响较小，1 mm填充层在剪切面附近的位移云线的面积较大，原因是填充层厚度小，填充颗粒在膨胀土颗粒挤压下发生错动较大，位移角大；在剪切第二阶段，颗粒的位移进一步增大，剪切面附近等值线开始密集，1 mm填充层位移等值线变化较为平缓，1.5 mm和2 mm填充层在靠近侧面墙体的部分等值线出现波动较大，这是由于剪切带厚度的增加，膨胀土颗粒的嵌入填充层运动的深度不变，因此出现一定的波动；在剪切的第三阶段，不同的剪切层厚度的在剪切面附近均出现绿色条状位移带，厚度随着填充层的厚度增加而增加，等值线在1 mm填充层时变化最多，在剪切面上部整体呈现圆弧状，下部出现有“鱼刺”状波动.

3.3 颗粒转角分析

图3 颗粒位移示意图Fig.3 Particle displacement diagram

颗粒转动的角度是反映颗粒运动的重要参数，以剪切位移0 mm为初始状态，分别统计颗粒各个阶段的旋转量，同时以设定的特定转动量（本次模拟试验设为50°）为界限，将转动量大于特定值的颗粒分别设定为红色，对于转动量小于特定值的颗粒，根据其具体数值的大小再进行分类，结果如图4所示.可以发现，当填充层厚度为1 mm时，试样在剪切第一阶段，颗粒的转动主要集中在剪切面的附近，旋转的角度较小，并未形成连续贯通的剪切带；在剪切第二阶段，颗粒的旋转角度明显变大，剪切带开始贯通，周围土体颗粒的旋转角度开始增加，旋转的颗粒数量开始明显增加，旋转区域主要集中在剪切盒的中间位置，在剪切盒的左上角和右下角颗粒的转动不明显；在剪切第三阶段，剪切面附近带状区域开始完全贯通，剪切面周围土体形成较多带状区域.填充层厚度为1.5 mm和2 mm颗粒转动规律较为接近，在剪切第一阶段，颗粒就形成贯通的剪切带，随着剪切的进行颗粒的转动角度不断增加，剪切面附近剪切带的厚度随着填充层厚度的增加而明显增加，这是由于填充层颗粒的直径小于周围膨胀土颗粒的直径，所以颗粒转角在剪切面附近的转角最大；在剪切第二阶段，剪切面周围条状剪切带宽度明显增加，转角开始逐渐大于50°，剪切盒上部和下部块状区域开始贯通，剪切盒顶部颗粒的形成近似水平的剪切带；在剪切的第三阶段，剪切面附近颗粒的转角基本开始大于50°，呈现红色，但红色区域的面积近似相等，这表明填充层厚度大于1.5 mm后，剪切面周围的性质开始由剪切面附近的颗粒属性控制，上部颗粒和下部颗粒很难越过填充层颗粒.

为了更好地展示颗粒在剪切盒转动角度的变化，统计剪切带不同位置上下各0.5 mm范围内颗粒转角，求出旋转角度绝对值的平均值，得到颗粒旋转角度随试样位置和填充层厚度的变化规律，如图5所示.可以发现，当填充层厚度一定时，随着剪切的进行，颗粒在不同位置平均转角都在不断增大，尤其是在8～12 mm之间，颗粒在10 mm附近颗粒的平均转角达到最大值，剪切位移相同的情况下，颗粒的平均转角最大值随着填充层厚度的增加而减小，1 mm填充层的平均转角最大达到250°.而在高度为11 mm和9 mm附近时，颗粒的平均转角随着填充层厚度的增加而明显增加，原因是填充层颗粒的直径小于膨胀土颗粒，转动角度较大.

图4 颗粒转角变化示意图Fig.4 Schematic diagram of particle rotation angle change

3.4 力链

力链是颗粒体系内部接触互相连接逐渐形成能够传递荷载的通道，这些通道错综复杂地交织在颗粒体系内部形成力链网络［8］.在Pfc颗粒流程序中直接提供了力链选项卡，可以直接清楚地显示颗粒之间的力链变化，本文选取1 mm填充层和2 mm填充层的颗粒力链网络来研究土体细观的土骨架的受力情况.图6为试样力链变化规律，其中力链的宽度代表了颗粒之间作用力的大小.可以发现，土体在初始状态时，力链的方向基本分布在竖直方向，填充层由于颗粒数量较多，力链较为密集，在上部和下部强力链和弱力链交叉分布.随着剪切的进行，力链的方向逐渐开始沿逆时针方向倾斜.在剪切第二阶段，倾角与水平方向夹角接近45°，强力链出现明显的平行分布；在剪切第三阶段，强力链的倾角比第二阶段小，局部的强力链消失.2 mm填充层的力链数量明显多于1 mm填充层，且力链的数值也大于1 mm填充层数值.

4 结论

本文通过Pfc2D颗粒流软件对膨胀土裂隙填充的直剪试验进行了模拟，分析宏观和细观结果，主要结论如下：

1）数值试验和室内试验的结果较为接近，二者应力应变关系变化的趋势基本一致，模拟效果较好；

2）随着剪切的进行，颗粒的位移在不断增大，不同填充层厚度的颗粒位移剪切面附近的大小一致，但在剪切面上部和下部有明显差别，1 mm填充层的位移云图比其他情况变化均匀；

3）在剪切开始时剪切面的上部和下部出现局部不连续的剪切区域，随着剪切位移的增大，局部的旋转区域开始连接形成旋转带，填充层厚度为1.5 mm和2 mm的颗粒转动规律较为接近，1 mm填充层的转角相比其他情况最小；

4）颗粒的转角在剪切面附近达到峰值，最大达到250°，沿着高度向上部和下部逐渐减小，在7 mm以下和13 mm以上区域平均转角变化较小；

5）力链在剪切过程中变化明显，力链随着剪切的进行开始沿着逆时针方向发生旋转，在剪切应力达到峰值时，力链方向与水平方向的夹角约为45°，峰值后强力链在局部变化较大，剪切面附近的弱力链由于颗粒的数量较多，因此变化不明显.