边坡土体变形和剪切带形成机制室内试验与数值模拟研究

2021-08-25刘洪波李水江叶永康梁崇旭

广东土木与建筑 2021年8期

刘洪波，李水江，叶永康，梁崇旭

（1、广州环投增城环保能源有限公司广州511335；2、广州环保投资集团有限公司广州510330；3、上海大学土木工程系上海200444）

0 引言

我国是一个地质环境脆弱、地质灾害频发的国家，而滑坡是众多地质灾害中发生最为频繁的类型，2018年全国共发生地质灾害2 966起，滑坡就占据了1 631起［1］。

滑波是坡体因多种因素耦合变形，最终被某些诱发因素（土体抗剪强度减小、降雨入渗、土压力增大等）激发失稳、产生滑动的一种地质灾害现象。目前在工程上对滑坡的稳定性分析主要采用极限平衡法计算稳定系数进行判定，但极限平衡条分法在计算原理上对滑坡体进行了诸多简化，将滑坡岩土体假定为刚性，难以反映滑体内部应力应变特征，并且对滑坡剪切带的轨迹进行人为的规定，难以反映滑动面（剪切带）触发及扩展机制，边坡渐进性破坏过程难以深入研究［2］。滑坡计算也可采用数值模拟方法，其主要优势在于，能够分析滑坡变形破坏过程中的滑体运动特征及应力－应变分布特征，确定滑坡变形破坏机制，并可以求解较为符合真实情况的稳定系数，提高效率。

WOLF等人［3］在拉伸荷载和离心机增加重力作用下，使用砂土进行了离心模型试验，研究了不同应力水平条件对剪切带形成、剪切带厚度以及破坏面与最小应力夹角等的影响。陈强等人［4］通过离心试验，研究了多级平台砂土边坡在不同含水率下的边坡稳定性差异，并对边坡破坏的机理进行总结等。刘君等人［5］探索了数字图像技术在大型振动台模型试验中应用的可能性，笔者等人使用细粉砂制作边坡模型，输入频率10 Hz的水平地震波，采集试验全过程的边坡侧面图像。RECHENMACHER等人［6-7］采用平面应变试验、三轴试验和数字图像相关技术测量了试样的局部位移特征，量化了剪切带体应变、剪切带倾角和厚度等物理特性，并分析了不同颗粒直径的砂土试样得到的试验结果差异。孔亮等人［8］基于数字图像相关技术和直剪试验，通过改造直剪仪开发了一种砂土细观测量系统，使用图像处理软件计算试验过程采集的图像，得到了试验不同时刻任意测点的位移场、应变场与颗粒方向角等数据。薛杰等人［9］使用改进的三轴试验数字图像相关测量系统，选取福建标准砂在不同的围压条件下进行三轴试验，分析了不同工况的应力-应变关系，通过试验全过程的应变场的变化特征研究了剪切带的触发及扩展情况。

现有文献很多集中于数值模拟研究滑坡内部剪切带、变形破坏和稳定性的变化。宋浩燃等人［10］以某隧道弃渣边坡为研究对象，利用三维离散元软件PFC3D对弃渣边坡变形至滑动过程进行数值模拟，研究其变形过程，结果表明边坡的变形破坏主要集中在坡体表层，滑动的规模受弃渣边坡角度影响；李巧学等人［11］运用三维离散元程序，对已失稳边坡还原地质体并进行数值计算分析，阐明岩质边坡在失稳破坏过程中速度和位移的变化规律；徐光霁等人［12］采用离散单元法，研究黏性土边坡的安全系数和失稳过程；黄达等人［13］进一步探明软硬互层反倾边坡的倾倒变形机制，融合离心模型试验与UDEC模拟，研究此类边坡的破坏模式与影响因素，并通过点对分析，讨论了变形的力学机制，在岩层内预置随机裂隙，获得了破裂面的演化规律；曾红艳等人［14］利用GeoStudio有限元软件，探究降雨入渗引起的膨胀土边坡变形破坏机理和在不同竖向植筋带长度、设置间距等影响因素下膨胀土边坡的变形破坏特征；陈观明［15］通过模拟分析顺层岩质边坡应力场、位移场及剪切区分布情况，研究了顺层岩质边坡的滑坡机理；毛广志等人［16］基于FLAC3D开展路堑边坡楔形体破坏机理数值模拟研究；杨文琦等人［17］在传统剩余滑坡推力法中引入膨胀力项，得出了适用于膨胀土边坡的滑坡稳定性系数计算公式，并通过离心试验和模拟验证了该计算公式；左健扬等人［18］基于Petrel软件实现的三维地质模型，提出了一种考虑整体滑动的三维极限平衡新方法，并通过经典算例的验证，证实该法的可行性和实用性；胡谢飞等人［19］将岩土体物理力学参数取值区间数学与有限元方法结合起来，采用区间有限元法对边坡稳定性进行不确定分析，编制区间有限元程序，结果表明重度、粘聚力及内摩擦角区间对安全系数区间的影响较大；HU等人［20］通过在重力增加过程中向摩擦角添加校正因子，开发了可以消除GIM中误差的改进重力增加方法（MGIM）。通过将GIM，剪切强度折减法（SSR）和MGIM应用于不同尺寸下的均质边坡模型，对MGIM进行了比较研究，结果表明MGIM能较好地反映土坡土工离心模型试验的变形破坏特征；CHEN等人［21］提出了一种基于欧拉-拉格朗日耦合（CEL）数值策略和基于能量的准则的边坡稳定性分析方法，对于软带较薄的边坡，本方法可以获得局部破坏和整体破坏的可靠结果；朱大勇等人［22］根据最优控制理论提出了边坡临界滑动场的概念，并提出了模拟临界滑动场的数值方法，算例表明临界滑动场中的临界滑动面是理论解的数值逼近。目前研究成果主要集中在边坡的数值模拟，未能和实际试验对照研究，试验与模拟未能有机统一，缺乏一定的实际意义，对于边坡的失稳破坏只起到辅助说明的作用，

本文采用DIC图像分析技术与室内模型试验结合，分析边坡土体变形和剪切带形成机制。同时基于有限差分软件FLAC3D，对边坡渐进性破坏全过程下的内部剪切带、变形破坏机制以及稳定性的变化进行研究。最后结合室内试验和数值模拟结果对边坡失稳破坏的机理进行研究。

1 室内试验和DIC分析

1.1 室内试验模型

砂土边坡模型是由福建标准砂制备，按照5%的含水率将砂土和水用搅拌机搅拌均匀，采用台阶式铺填方法，雨落法分层刮毛填筑，按照设计位置预埋传感器，如图1所示，逐层夯实保持模型密实度均匀。

图1 模型尺寸示意图Fig.1 Schematic of Test Model（mm）

模型试验使用万能材料试验机CTM进行加载，加载速率为4 N/s，并通过条形加载板进行力的加载，加载板长为纵向长为500 mm，横截面宽为100 mm，厚度40 mm。

模型尺寸为600 mm×350 mm×290 mm，坡脚角度为40°，边坡下方基层厚度100 mm，边坡高度250 mm，坡面长度为200 mm。

1.2 DIC测量系统

DIC图像采集系统主要由图像采集设备、光源和计算机组成，如图2所示。图像采集设备常用CCD相机或者单反相机，负责采集试验过程中的视频或者图片。摄像补光灯给试样表面提供均匀的光照。计算机编程处理采集的图像，进行相关计算以得到试样的变形信息。本研究采用自主编制的MATLAB程序进行相关性计算。

图2 DIC测量系统Fig.2 Measure System of DIC

1.3 坡顶荷载-沉降曲线

试验使用万能材料试验机进行加载。加载速度设定为4 N/s，直至砂土边坡模型破坏。万能材料试验机可以自动记录试验过程中坡顶荷载P、坡顶位移数据，位移以竖直向下为正。40°坡角砂土边坡P-s曲线如图3所示。

图3 40°坡角砂土边坡P-s曲线Fig.3 P-s Curve of 40°Sand Slope

P-s曲线变化分为2个阶段，第一阶段，伴随荷载增大，坡顶沉降接近线性增大，坡顶砂土处于整体压密状态，还未产生剪切破坏。第二阶段，沉降继续增大，竖向荷载反而下降，此时砂土边坡发生破坏失稳，剪切带内土体颗粒会沿着滑动面开始向下滑动，边坡失去承载能力。40°坡角砂土边坡极限承载力约为217 kPa。

1.4 边坡位移云图分析

将试验过程不同时刻的照片与加载初始时刻的照片进行MATLAB计算。T是加载初始时刻至边坡失稳破坏（边坡达到极限承载力）的总时间，水平位移以向右为正，竖向位移以向下为正。

由图4可知，加载初期，水平位移较小，主要出现在加载板右侧土体的浅层区域。随着加载持续，颗粒之间位置调整、应力向下传递，土体受影响的范围向右下方扩大，加载板右侧的土体向临坡面运动；加载板左侧出现较小向左的水平位移。破坏时刻，水平位移最大值出现在上部区域，并且由坡顶至坡底，水平位移逐渐减小，与模拟结果相类似。结合图3竖向荷载P-s沉降曲线），砂土边坡的破坏具有较明显的突发性，边坡出现由坡顶贯通至坡趾上方的的滑动面，滑动面上部接近圆弧状，下部水平呈直线型。

图4 基于DIC技术40°坡角砂土边坡水平、竖向位移云图Fig.4 Horizontal and Vertical Displacement Diagram of 40°Sand Slope during Loading Failure Based on DIC Technology（mm）

如图4所示，加载初期，竖向位移主要出现在加载板下方，说明竖向荷载较小时，承压板下方土体表现为竖向的压密过程。随着加载持续，承压板下方土体竖向位移增大且向下发展，应力通过土体颗粒间向下传递。加载板中轴线以下竖向位移从上至下逐渐减小，应力传递逐渐减弱，影响范围逐渐扩大。由竖向位移云图4⒟可知，竖直方向受影响的土体范围较大，接近于模型箱的边界。

1.5 边坡矢量位移云图分析

如图5所示，矢量位移图可以有效地观察到边坡受载以后变形、破坏的过程。加载板下方土体产生沉降，并向两侧产生斜向下的运动。因为左侧模型箱左边界受限的原因，土体主要向右侧临坡面运动，模拟左侧土体位移太小，位移场走势显示不是很明显，可能是受到左侧边界约束条件的影响；加载板右侧较浅的土层产生较小的向上的竖直位移，引起一定的坡面隆起，随着加载持续，竖向荷载越来越大，应力通过粒链向下传递，但传递逐渐减弱，土体受影响的范围持续增大。加载板右侧浅层土体产生角度（与水平线的夹角）较大的斜向右下方的位移，土体位置越深，角度越小，而临坡面的土体是以较小的角度向右上方运动。最终在破坏时刻，形成上部圆弧，下部接近水平的滑动面，边坡完全失稳破坏。

图5 基于DIC技术40°坡角砂土边坡加载破坏矢量位移Fig.5 Vetor Displacement Diagram of 40°Sand Slope during Loading Failure Based on DIC Technology

1.6 边坡剪应变云图分析

边坡在达到极限承载力后会失稳发生破坏，临坡面土体会沿着某一断面（即滑动面）突然滑动，这种断面称之为剪切带，而剪切带的形成往往伴随着应变局部化现象的出现。40°坡角砂土边坡加载破坏时刻剪应变云图如图6所示，横纵坐标为模型实际尺寸。从剪应变云图可以更直观的观察到剪切带贯通形成圆弧滑动面的过程。

由图6可知，加载板两侧及坡顶浅层区域出现较小的剪应变。随着加载的持续，出现不均匀应变的区域向土体深处延伸，加载板两侧出现了明显的剪切变形，土层上部加载板两侧的剪应变有连通在一起形成连续剪切面的趋势，此时土体内正处于剪切带的初期形成阶段。破坏时刻，DIC试验中加载板两侧形成2条剪切带在坡中位置连接在一起，形成1条从坡顶贯通至坡脚的剪切带。同时，在剪切带上剪应变基本为正值，剪应变大小差异较大，分析其原因，可能是剪切带上砂颗粒之间摩擦发生位置重构，砂颗粒间孔隙变大，局部位置砂土颗粒发生明显的扩张，剪应变较大。

图6 基于DIC技术40°坡角砂土边坡加载破坏剪应变云图Fig.6 Strain Diagram of 40°Sand Slope during Loading Failure Based on Dic technology

坡体上部局部位置剪切带倾角约为115°，接近坡体下部的局部位置二剪切带倾角约为150°，剪切带与水平线的夹角从上至下逐渐增大。并且发现本次试验40°坡角砂土边坡剪切带厚度约为3～5 mm（6～10倍平均粒径）。

1.7 边坡土体土压力分析

由图7可知，6个测点土压力变化曲线具有相似的规律。在土体破坏失稳前经历缓慢增大、快速增大2个阶段。在土体失稳以后，各监测点土压力急速下降至某一稳定值。2#、5#、8#监测点位于加载板中线位置，该位置土压力从坡顶至坡底随着深度的增加越来越小。4#、5#、6#监测点位于边坡中部，5#监测点在加载板中线位置，此处土压力最大，两侧土压力较小，其中，在土压力值稳定以后更接近临坡面的6#测点数值最小。这是因为此处土体变形较大，土压力发生了卸载。

图7 试验中40°坡角砂土边坡不同测点土压力时程曲线Fig.7 Pressure Time History Curve of 40°Sand Slope in Tests

2 有限元数值模拟

2.1 FLAC3D模型

基于有限差分软件FLAC3D建立砂土边坡三维模型。共创建了25 875个实体单元，28 704个节点，计算模型示意图如图8所示，模型相关参数如下：砂土层，密度为1.74 g/cm3，重度为18.0 N/m3，体积模量为20 MPa，粘聚力为1 kPa，内摩擦角为31.2°。

图8 边坡计算模型Fig.8 Calculation Model of Slope

模型底部边界固定X、Y、Z方向速度约束；前后边界固定X方向速度约束；两侧边界固定Y方向速度约束；地表自由。

2.2 边坡位移云图分析

根据图9⒜超荷载法加载边坡破坏时的水平位移云图，可以看出水平位移最大值区域集中在上部临坡面，边坡表面发生一定的隆起，并且与加载板之间还存在一定厚度的土体，可能是由于与加载板接触的土体由于摩擦力和受荷密实的约束，从而将应力传递给与加载板间接接触的下部土体，从而水平位移最大值区域发生在受荷面下部一定距离处，并且由坡顶至坡底，水平位移逐渐减小。

根据图9⒝超荷载法加载边坡破坏时的竖向位移云图，可以看出边坡破坏时，靠近加载板下方的土体竖向位移最大，并不断向右下方发展，通过应力链将应力进行扩散，由于坡面没有边界约束，发生应力卸载，边坡中上表面局部发生一定向上的竖向位移，从而表现出一定的隆起现象。

图9 基于FLAC3D数值模拟40°坡角砂土边坡水平、竖向位移云图Fig.9 Horizontal and Vertical Displacement Diagram of 40°Sand Slope during Loading Failure in FLAC3D

2.3 边坡矢量位移云图分析

如图10所示，整体矢量位移的趋势与室内试验结果相似（见图5）。位移最大值区域集中在加载板下方，并向两侧产生斜向下的运动，有所不同的是模拟中左侧边界土体的位移较小，可能是受到左侧边界约束的结果，边界效应较大。位移场走势延深度方向和右侧临坡面方向，与竖直方向呈较小的夹角，坡面土体产生一定隆起，但不明显，这是由于土体粘聚力的影响。根据图5与图10的对比分析，室内模型箱试验和数值模拟得到的结果类似。

图10 40°坡角砂土边坡加载破坏矢量位移Fig.10 Vector Displacement of 40°Sand Slope during Loading Failure

2.4 边坡剪应变云图分析

如图11所示，边坡模型整体剪应变区域和室内试验结果（见图6）相似，两者的剪应变都是开始于加载板两侧及坡顶面浅层区域，随着不断加载，剪应变不断向深度方向延伸，且无论是FLAC3D模拟还是DIC室内试验，土层上部加载板两侧的剪应变都有连通在一起形成连续剪切面的趋势，FLAC3D模拟的剪切带主要由左侧加载板边缘产生，右侧边缘的剪切变形影响较小，基于DIC技术的室内试验有着相似的结果，两条从加载板两侧产生的剪切带在边坡中部连接，逐渐形成一条从坡顶贯通至坡脚的剪切带。临坡面通过应力卸载，在剪切面处产生较小的整体滑动，从而使滑动面土体颗粒间错动加剧，产生较大的正向剪切应变。

图11 40°坡角砂土边坡加载破坏剪应变云图Fig.11 Strain Diagram of 40°Sand Slope during Loading Failure

2.5 边坡土体内部应力分析

由图12所示，无论是水平还是竖向应力，加载板附近始终是土应力最大值区域；除了加载板下方，临坡面和左侧边界处土体应力比较小，因为土体受到挤压，土体颗粒通过应力链向下方和两侧传递，随着深度和水平距离增加，应力传递逐渐减弱，接近于临坡面的测点因为应力卸载的原因，土应力最小；而竖向应力变化主要集中在加载板中轴线附近，随着深度增加，应力逐渐减小，影响范围不断向坡底延伸。

图12 40°坡角砂土边坡加载破坏土中应力云图Fig.12 Stress Diagram of 40°Sand Slope during Loading Failure

由图13可知，模型内部土应力发展规律和室内试验结果相一致。随着不断加载，土中应力不断增大，且变化率逐渐增大，边坡破坏后，各监测点土应力骤降至稳定。根据2#测点和8#测点、4#测点和6#测点应力值范围，说明加载板轴线处应力延深度和两侧逐渐减小，边坡中部水平线处应力，以轴线处5#测点为原点向两侧减小，应力传递至临坡面，应力值达到最小值，这与室内试验结果（见图8）和模拟的应力云图规律一致。两者只在部分位置存在差异，例如边坡最终破坏时的土压力都比室内试验要小，可能是试验加载系统使得土体压密，土中应力提高，且边坡破坏时刻附近时间区域中，在同一水平线上的土压力数值基本相同，而数值模拟的结果存在一定的不同，这是由于试验数据采集系统的误差导致的。

图13 试验和模拟40°坡角砂土边坡土压力时程曲线Fig.13 Test and Simulation Time History Curve of Soil Pressure in 40°Sand Slope during Loading Failure

2.6 边坡三维破坏分析

将砂土边坡模型在加载破坏时刻时的位移场、位移矢量云图、应力场和剪切应变云图进行切片处理，能够直观的观察到边坡模型Y方向（纵向）的变形破坏状态。

由图14可知，无论是位移场，应力场变化还是剪应变区域，剪切带产生位置，在Y方向的数据结果都十分相似，产生的差异较小。土体变形主要集中在加载板轴线附近并向深度和两侧扩展。应力延深度和水平方向扩展，逐渐削弱，临坡面应力卸载，应力值最小。剪切带产生于加载板两侧边缘，左侧较为明显，最终由左侧剪切带主导形成贯通至坡脚的滑动面。