基于蠕变模型的边坡抗滑桩体系变形时效规律数值研究

黄凯湘

(广州市城市规划勘测设计研究院，广东 广州 510060)

1 概 述

边坡工程是工程地质十分重要的研究方向，同时也是环境地质、地质灾害、岩土工程的重要研究内容。随着国民经济飞速发展，人类对自然环境的改造、利用日益广泛，城市建设、铁路、公路及水利枢纽的修建，露天开采矿产资源、煤炭资源等，都迫切需要进行边坡工程的研究和整治[1]，抗滑桩是治理边坡最有效的手段之一。许多边坡的破坏、抗滑桩抗滑效果、挡土墙墙后土压力的变化以及挡土墙的位移等现象都显示出边坡土体具有明显的流变特性[2]。目前，在传统的边坡稳定性分析和抗滑桩设计中，通常采用极限平衡理论，把岩土体材料视为刚弹塑性体，力学性质方面仅考虑其弹塑性质，釆用滑动面上的抗滑力和下滑力之比来评价边坡的安全稳定性和进行支挡设计，这种稳定分析方法忽略了岩土体的蠕变特性，不能考虑边坡随时间发展，蠕变变形增加，工程地质情况变化以及物理力学参数弱化，抗剪强度降低等问题，从而导致的边坡稳定性系数降低甚至整体滑动破坏。因此，了解土体的蠕变特性，并分析边坡抗滑桩体系随时间发展的变形破坏规律很有必要[3～9]。

2 数值模型建立

2.1 理想几何模型

如图1所示，有一无限长的土质边坡采用抗滑桩加固，坡高 10.0 m，坡度为1∶1.5，桩位置距离坡脚为 10.5 m，桩长 15.5 m，桩直径 0.8 m，桩间距为4倍桩径 3.2 m，桩端距离土体底部 2.0 m。无桩时，采用边坡剖面进行分析；有抗滑桩存在时，本问题不能简化为平面应变问题分析，这里利用对称性，取图中的阴影部分进行分析。设置5个监测点，分别设置在坡肩后缘、坡肩、坡体临空面中间、坡脚、坡脚前缘。

图1 边坡抗滑桩模型示意图

2.2 网格剖分

整个结构体系均采用实体单元，主要用到的单元有：C3D8(8节点六面体积分单元)，网格划分技术为sweep，单个网格直径为 0.4 m。土体计算模型共划分单元 27 252个，每根桩划分单元320个，如图2所示。

图2 边坡抗滑桩体系网格剖分图

2.3 土体强度和蠕变参数选取

(1)蠕变本构模型Singh-Mitchell模型介绍

分析中采用摩尔库伦弹塑性模型和Drucker-Prager与Singh-Mitchell弹塑蠕变耦合模型。从土流变的宏观角度出发建立的流变本构模型主要包括：理论模型、经验蠕变模型以及理论与实践相结合的半经验理论模型。理论模型无法描述软黏土的非线性蠕变，而经验蠕变模型的通用性不强，为了能更好拟合实际土体的蠕变特性，于是提出了一种依据理论模型的本构关系，并根据土体的试验结果而抽象的半经验理论模型。其中常用的半经验理论模型包括Mesri蠕变模型和Singh-Mitchell蠕变模型等，ABAQUS软件也提供了Singh-Mitchell蠕变模型[10]。因此，在本文的桩-土相互作用的蠕变模型中，土体的蠕变特性将通过Singh-Mitchell蠕变模型来体现。

(1)

(2)

(3)

则：

(4)

这样用方程确定应变，需要4个参数ɛ0，m，α，αA。

如果令ɛ0=0，方程变成

(5)

由于式(4)及式(5)均为描述土体蠕变的一种近似模型，而式(5)具有更简单的形式，通常把式(5)称为Singh-Mitchell三参数应力-应变-时间关系模型，在国内岩土工程中应用比较广泛。

(2)参数选取

本文采取某滑坡滑体碎石土土样，通过土体的常规土工试验、常规三轴及三轴流变试验，得出土的强度参数和蠕变特征，并利用Singh-Mitchell模型较好地拟合其在不同应力状态下的蠕变曲线，得到蠕变参数，如表1～表3所示。

边坡土体摩尔库伦模型参数 表1

边坡土体Druker-prager模型参数 表2

边坡土体Singh-Mitchel蠕变模型参数 表3

2.4 接触面模型

本文采用的接触模拟方法为罚函数法。接触面定义采用主-从接触算法，主控面上的节点可以穿透从属面，但从属面上的节点则不能穿透主控面，取桩表面为主控面，桩侧土体表面为从属面。接触面间的相互作用包括法向和切向，本模型桩-土接触面采用的是硬接触，即两物体只有在压紧状态时才能传递法向压力，若两物体之间有间隙时不传递法向压力。当接触面处于闭合状态时，接触面之间可以传递切向应力，即摩擦力，该模型的方程可表示为：

τcrit=μp

式中p表示法向接触压力；μ是摩擦系数，可以随剪切速率或其他场变量变化；

τcrit是摩擦力某一极限值，ABAQUS认为当摩擦力小于这一极限值时，接触面处于黏结状态，若摩擦力大于改极限值之后，接触面开始出现相对滑动变形，称滑移状态。本论文的桩土之间的摩擦系数取值为μ=tanφ=0.4。

2.5 边界条件

x、y、z方向详见图3、图4。限定模型左右两面上x向的位移，限定前后两面上y向的位移和底部x、y、z三个方向的位移，在桩的对称面上也要约束y方向上的位移。本模型为截取边坡抗滑桩体系中的一个桩土对称单元进行分析，对称面上的垂直轴向(y方向)位移为0，边坡底部和边坡前后缘边缘远离滑动变形区，约束相应方向变形，符合实际情况。

2.6 计算分析步

第一步，对土体施加的重力将与初始地应力正好平衡，使得土体的初始位移接近零；第二步，进行弹塑性变形受力特性分析；第三步，进行蠕变分析，蠕变模型的模拟时间为50年(43 800 h)。

3 边坡抗滑桩体系变形时效特征数值分析结果

3.1 蠕变效应下塑性区分布

经逐小时计算，最长时间达50年共 43 800 h，常规强度下的弹塑性模型、蠕变1年、30年、蠕变50年后的塑性区云图如图3所示。弹塑性模型下，只在抗滑桩底部出现微量的塑性区，边坡稳定；当考虑蠕变作用，塑性区首先出现在坡脚并随着时间推移往坡肩发展，在30年后边坡抗滑桩体系出现塑性贯通区。这是因为在坡脚的土体主要受上覆土体的重力“推挤”作用，随着时间的演变，土体流变效应逐步发挥作用，使土体本身的强度参数逐步弱化，变形逐渐增大，产生较小的剪切塑性破坏区。以塑性应变从坡脚到坡顶是否贯通作为边坡破坏的判据，则该加桩边坡在蠕变30年后发生破坏。

图3 不同时刻加桩边坡塑性区云图

3.2 蠕变效应下位移分布

如图4所示是蠕变50年时的位移等值线云图，抗滑桩阻止了其上部土体向下的滑动变形，没有出现整体圆弧状滑动面，但桩前土体仍然产生了向下的失稳滑动变形，该部分土体和桩是脱开的，桩后土体产生了绕桩滑动，即绕流失稳现象。

由图5不同时间坡表水平位移曲线可看出，初始弹塑性模型的位移大小相对于蠕变导致的位移大小可以忽略，蠕变1年后，由于蠕变初期位移增长较快，蠕变位移增量主要发生在x=8～x=28之间，即坡脚前 2 m和坡肩后 3 m之间的临空面，斜坡表层位移最大处增值约每20年增加 5 cm，即 0.25 cm/a，说明蠕变可能导致边坡临空面附近形成整体变形。

图4 加桩边坡蠕变50年位移云图(U单位/m)

图5 不同时间坡表水平位移曲线

图6 坡表监测点随时间x位移曲线

由图6坡表监测点随时间x位移曲线可看出：在流变初始阶段，边坡土体处于弹性变形阶段，位移量呈线性增加趋势，但这种增加趋势持续时间较短，体现为初始流变，然后，边坡土体逐渐进入到衰减蠕变阶段，蠕变位移增速逐渐减小直至为零，流变曲线逐渐趋于缓和，这表明边坡土体在经历一定时间的流变后，边坡逐渐进入到稳定流变状态。假如没有抗滑桩的存在，边坡土体在重力等作用下，坡肩位移应大于坡脚位移，本模型中由于抗滑桩的阻滑作用，桩后坡肩附近的1、2号监测点位移明显小于桩前的3、4号监测点，说明抗滑桩阻止了上部土体向下的滑动变形。

图7 坡肩2号监测点竖直深度水平位移

图8 坡中3号监测点竖直深度水平位移

由图7、图8竖直监测孔水平位移图可看出：在蠕变作用下，桩后坡肩2号监测点水平位移主要在深度 4 m以上，桩前坡中3号监测点水平位移主要在深度 10 m以上。这主要是因为抗滑桩降低了桩后土体的变形影响深度，使得桩后坡肩2号监测点变形深度小于桩前的3号监测点。

3.3 抗滑桩时效变形分析

图9给出了蠕变10年、蠕变50年桩的变形，由图可见桩的上部产生了弯曲变形，而下部则呈较好的锚固桩，是典型的中长桩变形模式。为了印证桩周土体的变形模式，将桩顶附近法向脱开变形绘制于图10，桩前土体在桩顶以下范围约2倍桩径的范围内产生了脱开变形，在蠕变50年后最终开裂 1.2 cm，桩后土体与桩之间都处于闭合状态。

图9 蠕变10年、50年桩的变形(U1单位/m)

图10 蠕变30年、50年桩顶土体开裂

图11为流变50年后桩身土压力，由其可知土体蠕滑50年后，抗滑桩上的土压力比常规分析中所假设的桩前被动土压力，桩后主动土压力要复杂得多。桩前顶部以下约2D范围内由于土体的滑动失稳，土体与桩之间是脱开的，因而没有土压力(接触压力)。在较深的土体中，桩前土体的土压力要大于桩后的数值，这是和桩的弯曲变形密切相关的。

图11 桩后应力

3.4 桩土相互作用时效分析

图12表明，边坡抗滑桩附近变形有明显的土拱效应，在蠕变稳定期如图13所示，蠕变50年比蠕变30年的位移增量仅为 1 cm，抗滑桩对于阻止蠕变变形有很好的效果。

图12 蠕变50年与30年间的位移增量

图13 桩前两米土体水平位移

图14为桩后两米横向最大主应力，由于网格在桩位置的非完全对称性，导致应力拱非完全对称，由图可知桩后 2 m形成了应力拱，且随着蠕变时间的推移，桩端所承受荷载增强，拱顶土体最大应力也有所增大，土拱作用得到增强。

图14 桩后2 m最大主应力

图15为抗滑桩与桩前土体协同变形时效曲线，桩前土体在土蠕变性质影响下随着时间产生了水平位移，而抗滑桩整桩置于碎石土上也随着碎石土整体蠕变产生了向坡脚方向的水平位移，抗滑桩位移在蠕变50年后位移为 10 cm。在蠕变两年后，桩顶土体与桩在地表开始出现开裂约为 2 mm；在碎石土蠕变作用下，随着时间发展桩土开裂不断变大并往土体深部扩展，在蠕变50年后地表桩土开裂变形为 12 mm，开裂深度约为 2 m。

图15 抗滑桩与桩前土体协同变形时效曲线

3.5 加桩前后边坡变形时效特征对比

(1)塑性变形对比

对无桩边坡进行数值模拟，无桩边坡在蠕变下塑性区贯通时间为20年(如图16)，此时坡脚最大塑性值为1.079e-1，坡肩最大塑性值为2.694e-2；而加桩边坡蠕变下塑性区基本贯通时间为30年(图3)，坡脚最大塑性值为9.89e-2，坡肩最大塑性值为1.645e-2。由此可知，抗滑桩有效地阻止了抗滑桩塑性区的贯通，阻止了边坡的整体滑动。

图16 无桩边坡塑性变形

(2)位移分布对比

为了分析抗滑桩对边坡碎石土蠕变作用的阻滑效果，将无桩边坡和加桩边坡蠕变50年表面土体位移绘制在图17中，加桩边坡主要阻止了坡脚到坡肩的斜坡部分的位移，使得最大位移由 18 cm降为 13 cm。

图17 加桩与否边坡蠕变50年位移云图

4 结 论

(1)Singh-Mitchell模型可较好地拟合碎石土的蠕变本构关系，通过ABAQUS软件中的Drucker-Prager与Singh-Mitchell弹塑蠕变耦合模型能较好模拟边坡抗滑桩的蠕变时效规律。土体常规弹塑性本构模型仅能用于分析某一时刻边坡的稳定性，而采用弹塑蠕变耦合模型，能分析边坡随时间推移的时效变化规律。

(2)蠕变效应使边坡抗滑桩体系的变形有显著增加，塑性变形首先在坡脚出现，随着时间推移逐步向坡肩发展，最大塑性变形位于坡脚。坡体位移主要集中在斜面位置，边坡土体的蠕变性质表现出初始蠕变、减速蠕变和稳定蠕变三个阶段，坡面位移最大往土体深部位移逐渐减小。

(3)抗滑桩用于治理蠕变变形效果明显，不仅显著减少了桩后土体的水平变形值，也降低了桩后土体的变形影响深度。蠕变50年后，无桩边坡最大位移为 18 cm，加桩边坡最大位移降为 13 cm。无桩边坡的塑性区贯通时间为20年，加桩边坡塑性区的贯通时间延长至30年。

(4)抗滑桩桩前土体仍然产生了向下的失稳滑动变形，该部分土体和桩是脱开的，桩后土体产生了绕桩滑动，即绕流失稳现象。在碎石土蠕变作用2年后，桩前顶部以下约2倍桩径范围内由于土体的蠕滑失稳，土体与桩之间脱开约 2 mm，随着时间发展桩土脱开不断变大并往土体深部发展，在蠕变50年后地表桩土开裂变形为 12 mm，开裂深度约为 2 m。

(5)抗滑桩与边坡碎石土间有明显的土拱效应，在碎石土蠕变作用下抗滑桩所承受荷载有所增加，抗滑桩的上部产生了弯曲变形，而下部则呈较好的锚固状，是典型的中长桩变形模式。土体蠕滑50年后，抗滑桩上的土压力比常规分析中所假设的桩前被动土压力，桩后主动土压力要复杂得多，桩前顶部以下约 2 m范围内由于土体的滑动失稳，土体与桩之间是脱开的，因而没有土压力，在较深的土体中，桩前土体的土压力要大于桩后的数值，这是和桩的弯曲变形密切相关的。