基于颗粒流方法的含粗粒滑带土宏细观力学特性

李泽闯，刘志斌，程培峰，张 昊，蔡启源

（东北林业大学 土木工程学院，黑龙江 哈尔滨 150040）

0 引言

在滑坡形成及变形机理方面，滑坡的变形、失稳与滑带土应力状态的改变、强度的减小有着紧密联系，滑带土力学特性研究显得异常重要.国内外许多学者对滑带土宏观力学性质开展了较为深入的研究，目的是为了更好地理解和控制滑坡，也取得了丰富的成果.CHEN X P[1]等通过原位直剪、三轴和直剪试验研究了古滑坡滑带土的力学特性，指出滑带土应力-位移曲线呈现应变硬化，没有显著的峰值，剪切应力随着剪切位移的增加而逐渐减小.NIAN T K[2]等通过对重塑滑带土试样进行直剪试验研究，指出非饱和滑带土具有“软化”应力-应变特性，而饱和滑带土具有“硬化”曲线.LIAN B Q[3]等通过对滑带土进行环剪试验，指出不同含水率下的残余强度都近似与法向应力成正比，且随着含水率的增加，滑带土的残余粘聚力和摩擦角也随之减小.高旭[4]等指出不同应变速率下滑带土抗剪强度值差异明显.WANG L[5]等通过对滑带土进行环剪试验，指出滑带土残余强度随剪切速率的增加而降低，剪切位移和剪切速率的增加会改变剪切带的结构.武立清[6]等通过对砾质滑带土进行直剪试验，研究了砾石的质量分数、颗粒级配等因素对其力学特性的影响.ZOU Z X[7]等通过原位剪切试验，研究了黄土坡滑坡滑带土的剪切力学性质，指出滑带土中砾石的存在增加了滑带土的强度.李远耀[8]等通过对三峡库区滑带土的抗剪强度参数进行统计分析，指出该地区滑带土的内摩擦角随碎石质量含量的增加而增大.陈晓平[9]等通过对含粗粒滑带土的抗剪强度进行试验研究，指出滑带土的抗剪强度受到粗粒质量含量的影响.从文献[10]～文献[11]对滑带土的研究总结中可以看出，学者们逐渐发现大量滑坡滑带上的土一般都含有砾石等粗颗粒，含粗粒的滑带土（粗颗粒的质量分数为25 %～50 %）力学特性与一般细粒滑带土有较大不同，粗颗粒的质量占比是滑带土的变形和强度特性的重要影响因素.现阶段对于有较多含量（质量分数为25 %～50 %）粗颗粒的滑带土的研究较少，且关于滑带土力学性质的研究大都停留在宏观层面，宏观层面的研究不足以解释滑带土强度的控制因素及滑坡变形演化的破坏机理.鉴于此，一些学者借助扫描电镜、X 射线衍射仪等微观研究设备对滑带土的微观结构和矿物成分进行了深入的研究分析.刘动[12]等采用扫描电镜指出滑带土颗粒在剪切过程中有明显的定向排列现象，微观结构的变化与宏观力学特性有着密切联系.MIAO F S[13]等采用压汞试验研究了三峡库区滑带土的微观结构，指出微观的分形结构单元是影响滑带土力学性质的关键因素.但微观研究主要反映在以黏土矿物成分为主的细颗粒滑带土（大于0.075 mm 的颗粒质量分数小于25%）.受到粗颗粒尺寸的限制，微观试验的研究难以解释较多粗颗粒组分滑带土的特性，不适合对于含粗粒滑带土开展研究.因而对含粗粒滑带土进行更加精细和深入的研究出现了瓶颈.

目前，在滑带土方面的研究主要集中细粒滑带土的强度特性和微观特征，而针对含粗粒滑带土变形破坏机理等方面的研究还比较欠缺.可利用数值模拟手段突破颗粒尺寸的限制，进行含粗粒滑带土细观层面的研究，以便与滑带土宏观剪切试验研究相结合.滑带土的宏观力学参数实质上是土颗粒细观参数的外在表现，颗粒流方法是模拟土体细观力学的有效手段，可以通过颗粒流数值模拟研究含粗粒滑带土细观参数的变化对其宏观现象的影响.郑博宁[14]等提出对于含砾滑带土可以通过颗粒流方法及CT扫描技术建立三维模型.含粗粒滑带土的宏-细观研究尚处于起步阶段，目前研究成果较少，有关细观剪切变形特性的研究不够全面.因此，研究含粗粒滑带土的宏细观力学特性，既能从宏观层面理解和认识其力学特性，又可以在细观角度解释含粗粒滑带土的细观剪切变形机理，对于探究滑坡的变形破坏和稳定性分析具有重要意义.

采用可视化直剪和颗粒流数值模拟相结合的方法探究含粗粒滑带土的宏细观力学特性，先通过含粗粒滑带土的室内可视化直剪试验对其宏观力学行为进行分析，之后运用颗粒流数值模拟在细观层面上对含粗粒滑带土的颗粒接触力链、接触力分布、位移矢量场、内部应力分布等细观力学参量的演化过程进行对比分析.

1 含粗粒滑带土的室内试验

土体剪切试验所用仪器为可视化应变控制式直剪仪，见图1.

图1 可视化直剪试验仪Fig. 1 visual direct shear tester

以往土工试验中常用的是不透明钢板组成的封闭式剪切盒，后者仅可在直剪试验结束之后才能看到剪切的结果.与常规剪切盒相比，可视化剪切盒不仅可以直观的观察到土体内部剪切变化的过程，而且还可借助工业显微镜记录剪切的全过程，方便后续分析和处理试验数据.通过对所记录的剪切图片进行图像处理，可为建立含粗粒滑带土的细观参数与其宏观力学特性之间的力学响应提供帮助.

1.1 滑带土基本特征

以抚顺西露天矿南帮滑坡为研究对象，在滑坡后缘I 号地裂缝的断陷陡坎处采集含粗粒滑带土试样[15].通过试验测得取样点滑带土的天然含水质量分数为21%，密度为19.1 kN/m3，塑限为16.2%，液限为26.2%.通过颗粒分析试验获得其粒径级配曲线见图2.

图2 含粗粒滑带土级配Fig.2 grading of soil with coarse grained slip zone soil

图2中，粒径区间组成分别为：0.500～2.000 mm占比10%、0.250～0.500 mm 占比20%、0.075～0.250 mm 占比20%、小于0.075 mm 占比50%.

1.2 直剪试验及可视化分析

对试样分别进行100 kPa、200 kPa 和300 kPa三组法向压力下的可视化直剪试验，3 组法向压力下的剪应力-剪位移变化见图3，得到抗剪强度参数黏聚力为80.68 kPa，内摩擦角为24.04°.

图3 含粗粒滑带土室内直剪试验结果Fig.3 laboratory shear test results of coarse-grained slip zone soil

对比图3 不同法向压力下含粗粒滑带土的剪应力-剪位移可以看出，含粗粒滑带土的剪切应力随着法向压力的增大而增大，100 kPa 和200 kPa 下的剪应力差值略大于200 kPa 和300 kPa 下的差值，且三组不同法向压力下的含粗粒滑带土的剪应力-剪位移曲线均呈现出明显的应变硬化现象.

在进行含粗粒滑带土的直剪试验时，利用直剪盒的可视化，采用工业高清数码仪器记录下了剪切试验的全过程，直剪试验的记录结果见图4.

图4 直剪试验100 kPa 法向压力Fig. 4 direct shear test at 100 kPa normal press

对比图4（a）、图4（b），可在宏观层面直观地看出含粗粒滑带土中的粗颗粒在直剪前后有明显的颗粒迁移情况，即粗颗粒在剪切过程中会发生错动、翻转等一系列不规则运动.此外对比剪切前后可以发现含粗粒滑带土中剪切带的发展情况，剪切带从试样中的微小裂纹和粗颗粒附近开始，而后逐渐发展扩大和贯通，最终导致试样发生剪切破坏，可视化直剪在一定程度上也可以为岩土材料的剪切带和局部应变化理论的发展提供实验支撑.

2 含粗粒滑带土的颗粒流数值试验

2.1 颗粒流及成样方法

由于土体天然的非连续性和力学行为的复杂性，使用颗粒流方法研究土体更加贴近土体真实的物理力学性质，它通过数值方法将岩土材料划分为有代表性的数千乃至上万的颗粒单元，是分析岩土材料内部颗粒间接触力和应变的有效工具.此外颗粒流方法可以有效模拟介质的开裂、分离等非连续现象，进而在细观层面上揭示介质的损伤及破坏机理.许多学者运用颗粒流对边坡破坏、滑坡等工程问题开展系列研究，并取得很多研究成果[16-19]，但针对含粗粒滑带土力学特性的细观研究较少.笔者通过引入颗粒流对含粗粒滑带土的直剪试验进行数值模拟，并从位移矢量图、大小主应力分布变化、颗粒接触力链以及组构等细观角度剖析其力学性质.

图5为建立的颗粒流模型.该模型按照含粗粒滑带土真实颗粒级配曲线（见图2）创建，颗粒共33 641 个，可以在细观层面上真实反映出含粗粒滑带土的力学特性.颗粒流模拟过程分为成样、固结、加胶结、加载4 步进行.在生成试样之后利用颗粒流内嵌fish 语言编写伺服函数，保证模型上下墙体施加恒定的法向压力，而后根据试验要求，上剪切盒保持不动，通过赋予下剪切盒墙体以速度进行剪切，直至模型试样发生破坏.在施加载荷过程中，每剪切10%即0.2 mm 记录一次模型的位移矢量场、颗粒接触力链等参量的变化.

图5 含粗粒滑带土颗粒流模型Fig.5 particle flow model of soil with coarse grained slippage

颗粒流数值模拟中，颗粒的力学性质通过颗粒的接触参数、接触模型等细观参数来体现.当前阶段含粗粒滑带土的宏观物理力学参数还不能直接转化成颗粒的细观参数，通过“试错法”进行一系列的数值试验得到与含粗粒滑带土宏观力学参数吻合程度较好的1 组细观参数，见表1.

表1 颗粒流模型细观参数Tab.1 microscopic parameters of particle flow model

2.2 抗剪强度参数对比

考虑到所要研究的含粗粒滑带土具有一定的黏性，根据以往研究经验，在颗粒流建模时选取线性接触黏结模型（Linear Contact Bond Model），在一定程度上保证了模拟的合理性和精准性.室内试验和颗粒流数值试验的拟合结果见图6.

图6 3 种法向压力下颗粒流拟合结果Fig. 6 three sets of discrete element model fitting results under normal pressure

对比三组法向压力下的模拟结果可以发现，数值实验的剪应力-剪位移曲线的初始弹性模量略大于室内试验模型及曲线，但整体上模拟曲线和试验曲线趋势变化一致、拟合程度较好，表明所选取的颗粒流细观参数较为合理，可为后续含粗粒滑带土的颗粒流相关研究提供借鉴和参考.

将颗粒流数值试验获得的抗剪强度包络线与室内直剪试验的结果进行对比，见图7，并将二者的抗剪强度参数进行对比，列于表2.从图7 可以看出含粗粒滑带土的剪切峰值应力随着法向应力的增大而线性增大，且室内试验值和数值实验的抗剪强度包络线的趋势相同，拟合程度较好.从表2 可以看出模拟值与试验值的抗剪强度参数相差不大，表明本文颗粒流模拟的合理性，此外可以看出含粗粒滑带土属于黏性土，其抗剪强度受到内摩擦角和黏聚力的共同控制，且内摩擦角相比无黏性土而言较小.

图7 试验与模拟抗剪强度包线对比Fig. 7 comparison between test and simulated shear strength cladding

表2 试验与模拟抗剪强度参数对比Tab. 2 comparison of experimental and simulated shear strength parameters

2.3 含粗粒滑带土剪胀性分析

在颗粒流数值模拟过程中，通过颗粒流内嵌fish语言编写函数，监测并记录含粗粒滑带土颗粒流模型在100 kPa、200 kPa、300 kPa 下的法向位移随着剪切位移的变化的剪胀曲线，见图8.

图8 100～300 kPa 法向压力下数值试验剪胀曲线对比Fig.8 comparison of dilatation curves in numerical tests under normal pressure of 100～300 kPa

从图8 可以看出在各法向应力下含粗粒滑带土数值试样的剪胀曲线整体呈显著剪胀性，法向位移随法向压力的增大而减小；在相同剪切位移下，法向应力越小，对应的法向位移越大，剪胀性也更明显.此外，在剪切的初始阶段，对比3 组不同法向应力下剪胀曲线可以看出，法向应力处于较低法向应力（100 kPa、200 kPa）时，法向位移为正，呈现剪胀性；法向应力处于较高法向应力（300 kPa）时，法向位移为负，呈现剪缩性；之后随着剪切位移的增加，所有法向应力下的法向位移都为正值，呈现显著剪胀性，上述剪胀性规律与张强[20]等基于颗粒流以土石混合体为研究对象所得的剪胀性趋势一致.从细观层面分析含粗粒滑带土数值试样的剪胀性主要是由于剪切开始前数值试样粗细颗粒为均匀分布，在剪切过程中颗粒会产生相互挤压和碰撞，导致粗颗粒发生错动和翻转，进而引起法向位移的增大和数值试样的体积膨胀，表现出剪胀性.在低法向应力下，粗细颗粒更容易发生翻转和错动，数值试样体积膨胀较大，因此其对应的法向位移较大，剪胀性更显著，而在高法向应力下，粗细颗粒不易发生错动和翻转，故其法向位移随剪切位移的变化相对较小，剪胀性没有低法向应力下显著.

3 颗粒流数值试验细观力学分析

3.1 位移矢量场演化分析

以200 kPa 法向压力为例分析含粗粒滑带土各细观力学参量的演化过程.根据含粗粒滑带土模型的位移矢量场的分布演化，可以从细观层面分析其在剪切过程中剪切带的形成和演化情况.图9 依次是颗粒流模型在固结之后以及剪切位移达到20%、50%、80%和剪切完成时的位移矢量场分布.对比分析可以发现，在模型固结之后位移矢量场表现为无规则的随机分布，但在模型的粗颗粒附近的位移较大，而细小颗粒附近的位移则较小.之后随着剪切力的施加，含粗粒滑带土模型的位移场开始呈现出沿着剪切力的方向分布，剪切下盒的颗粒位移整体大于上盒并呈现出一定的规律性，即在剪切应变较小时（20 %），颗粒位移主要集中于模型的左下角和右下角，而左上角和右上角出现明显的位移区分带（裂隙），右侧区域更明显；之后随着剪位移的继续增加（50 %～80 %），颗粒位移普遍开始增大，模型上盒颗粒在下盒颗粒的挤压和带动之下也沿着剪切方向运动，但右上角颗粒位移呈现明显疏松状态. 颗粒位移从左下侧到右上侧逐渐减小，这主要是由于左下侧墙体为剪切力施加一侧，而右上侧墙体属于上盒，在剪切过程中保持不动所引起；而后随着剪切位移达到最大，剪切结束时颗粒位移场分带规律显著，从左下侧到右上侧数值从大到小呈曲边条带状分布，上述位移场演化过程与薛亚东[21]等以崩积混合体为研究对象，进行直剪试验PFC 数值模拟的研究成果相似.最后对比分析还可以发现随着剪位移的增加，模型的左上角和右下角分别向下和向上凸出，表现为明显的剪胀性.

图9 位移矢量场分布Fig. 9 displacement vector field distribution

3.2 颗粒接触力链演化分析

散体材料中，接触力经由颗粒间接触从一个颗粒传递至另一个颗粒，这一传力特征在可视化后表现为链状结构，直观地称为力链.力链分叉、合并、交叉形成复杂的接触力网，主导散体材料的力学性质[22].在众多细观特征参数中，力链作为特征颗粒组，在尺度上介于颗粒与单元体之间，同时也是颗粒材料承载和传力的主体，既表现出与宏观应力、应变相同的变化特征，又能从机制上解决宏观力学分析中出现的争议[23].颗粒接触力链能从细观角度展现含粗粒滑带土模型在直剪试验中颗粒受力、传力及其变化情况，力链分布及演化情况见图10.

图10 颗粒接触力链分布演化Fig.10 evolution of particle contact force chain distribution

从图10（a）可以发现，在试样固结之后其颗粒接触力链在各个方向上的空间分布较为随机，呈现非规律性，但模型粗颗粒附近的接触力链相较于细小颗粒而言力链较粗、数值较大（接触力链的粗细表示接触力的大小），可以认为含粗粒滑带土颗粒流模型刚开始是由大颗粒相互之间紧密接触构成土体骨架，从而起到传递和承担载荷的作用.而后随着下盒剪切力的施加可以发现，颗粒接触力链也开始沿着剪切力方向分布，并初步呈现一定的规律性（见图10（b）、图10（c）），即模型左上角和右下角区域的接触力链显著减少，变得稀疏，同时中间区域粗颗粒附近的力链也开始减小，而左下角和右上角分布的力链数和数值则显著增加，这主要是由于颗粒受到水平剪切力，数值试样的受力和传力骨架开始重分布、重构造.而后随着进一步的加载，上述颗粒接触力链分布规律更加明显（见图10（d）、图10（e）），颗粒接触力主要分布在以上盒右侧墙体和下盒左侧墙体所围成的四边形区域，但力链较之前增加较少，因为此时剪位移对应的剪应力已接近峰值剪应力，模型发生剪切破坏，上述接触力演化规律与杨升[24]等以土石混合体为研究对象，进行三维直剪的颗粒流模拟所得结论具有一致性.从上述细观分析中可以发现，土体试样细观接触力链的方向、粗细等变化对应着宏观层面上剪应力的变化情况，两者密切相关.

3.3 大小主应力演化

土体应力及其分布演化情况是反映土体细观力学的重要特性，为了更加清晰地了解含粗粒滑带土颗粒流模型在模拟过程中应力的变化，通过编写函数，以应力十字的方式记录及输出了不同剪切阶段下数值试样内部大小主应力的分布情况，见图11，分别是模型在剪切位移达到20%、50%、80%和剪切完成时的大小主应力分布，图中应力十字的长短表示应力的大小，其长轴和短轴分别代表最大和最小主应力.从图11（a）可以看出，在模型固结之后其应力的大小和方向整体分布比较均匀，但在局部即含粗粒滑带土的粗颗粒附近的应力值较大，而后随着剪切位移的增加（见图11（b）、图11（c）），模型内部应力分布迅速发生变化，在模型左上侧和右下侧区域的应力减小而左下侧和右上侧区域的应力增大，且中间区域粗颗粒附近的应力值也开始减小，这主要是由于剪切力的施加而引起的应力重分布.之后随着剪切位移的进一步增加可以看出（见图11（d）、图11（e）），应力分布规律更加明显，即在左下、右上、中间区域模型的应力都沿着土体所受剪切力的方向分布和传递，且在联通左下和右上两个区域的过渡区域，其应力分布为斜向上约30°的角度分布，上述过程在细观层面完整地展现出数值试样在剪切过程中内部大小主应力的分布及演化过程.

图11 大小主应力分布变化Fig. 11 variation of principal stress distribution

3.4 组构演化分析

为更加清晰地分析含粗粒滑带土在直剪过程中颗粒内部之间接触力的演化情况，对200 kPa 下不同剪切阶段（固结、剪位移达到20%、50%、80%、100%）的颗粒接触力进行统计分析，并分别绘制了颗粒法向和切向的接触力分布，见图12、图13.

图12 颗粒法向接触力统计分布Fig.12 statistical distribution diagram of tangential contact forces

图13 颗粒切向接触力统计分布Fig. 13 statistical distribution diagram of normal contact forces

由图12 可以发现，颗粒法向接触力在剪切刚开始时为圆形分布，之后随着剪切位移的增加，初始阶段的圆形分布开始发生变化，逐渐演化为沿30°与210°连线方向分布，变为椭圆形状分布，之后随着剪切位移的进一步增加，椭圆形分布的中间区域开始收缩变窄，而两端开始集中凸出，最终接触力呈现出花生形状分布.法向接触力的形状从圆形分布的各向同性，逐渐过渡到椭圆形分布，最后演化为花生状的各向异性分布.上述演化过程也更加直观地展现出含粗粒滑带土在直剪过程中由于在下剪切盒施加水平剪切力之后内部的法向接触力的细观演化情况.

由图13 可以发现，颗粒切向接触力在剪切过程的不同阶段大致沿着4 个方向分布，形如蝴蝶状.切向接触力在剪切的初始阶段左侧上下和右侧上下4个方向切向接触力的分布还较为分散，而后随着直剪位移的增加，4 个区域内的切向接触力各自逐渐合并连为一体，主要沿着60°与240°和150°与330°两条连线方向进行分布，且剪切位移越大，其4 个区域内的切向接触力的分布范围和面积越大.

4 结论

通过对含粗粒滑带土宏观层面的室内可视化直剪试验，以及细观层面颗粒流的数值模拟相结合的方法，探究了含粗粒滑带土在直剪过程中的变形破坏机理及其宏细观力学特性，得出的主要结论如下.

（1）含粗粒滑带土的剪应力-剪位移曲线整体呈应变硬化趋势，剪切过程中表现出明显的剪胀性.此外粗颗粒在剪切前后有明显的颗粒迁移变化，剪切带从试样中的微小裂纹和粗颗粒附近开始，而后逐渐发展扩大和贯通，最终导致试样发生剪切破坏.

（2）含粗粒滑带土的位移矢量场和接触力链在剪切开始时表现为无规则的随机分布，之后随剪位移的增加，位移场从左下侧到右上侧的数值逐渐减小，最终逐渐演化为曲边条带状分布，颗粒接触力链随剪位移的增加而沿着剪切力方向分布，粗力链骨架开始呈现，土体骨架被重新构造.

（3）在剪切的初始阶段，含粗粒滑带土内部大小主应力为均匀分布，之后在剪切力的作用下，内部应力发生重分布；法向接触力在初始剪切时为圆形分布，之后演化为椭圆形分布，最终呈现出花生状分布，切向接触力沿4 个方向分布，形如蝴蝶状，剪切位移越大，切向接触力分布面积越大.