膨胀各向异性对非饱和膨胀土边坡稳定性的影响

刘正楠，张锐，刘昭京，郑健龙

(1.长沙理工大学交通运输工程学院，湖南长沙，410114；2.中国铁路南宁局集团有限公司，广西南宁，530029)

膨胀土在全球分布十分广泛，因其富含蒙脱石矿物而表现出吸水膨胀变形的工程特性，且因土中片状黏土颗粒的高度定向而使其膨胀具有显著的各向异性[1-4]。当膨胀变形受到周围介质或构造物约束时，便会在土体内部产生附加膨胀力。监测结果表明，因侧向膨胀力的叠加作用，膨胀土的土压力将为上覆土自重应力和不计膨胀作用时的主动土压力的倍数[5-6]，导致膨胀土边坡发生局部被动破坏[7]，给公路、铁路、水利等行业的建设带来极大的危害和巨大的经济损失[8-10]。如何客观分析膨胀各向异性对边坡稳定性的影响，揭示降雨入渗条件下膨胀土增湿膨胀对边坡破坏的作用机理，对于膨胀土地区的工程建设具有较好的应用价值，同时也可为膨胀土支挡结构加固设计提供参考。

膨胀土边坡的失稳往往在雨季发生，降雨入渗作用将引起膨胀土边坡浅表层自坡脚开始出现渐进式破坏[11-12]，长期以来普遍认为这是基质吸力降低引起的土体抗剪强度下降所致[13-14]，干湿循环引起的浅表层开裂会加速这一过程。采用极限平衡法进行反算，其强度很低，有时甚至接近于土体的残余强度；干湿循环虽然可以引起土体强度降低，但并不足以导致膨胀土边坡失稳[15]。事实上，膨胀引起的应力重分布和变形也对膨胀土边坡的浅层破坏的发展起到非常重要的作用[16-17]。近年来，数值方法因具有可重复性、直观性等特性，已广泛用于膨胀土边坡稳定性分析。分析方法可大致分为3 种：第一种是先由饱和-非饱和渗流分析得到渗流场，再由边坡内部湿度分布规律结合本构模型求解[17]；第二种是通过建立流-固耦合方程将非饱和渗流场与应力场相互耦合求解[18]；第三种是利用热传导热能量平衡方程与孔隙渗流连续方程数学描述的相似性，推导出热传导膨胀模拟增湿膨胀的替代方程，结合本构方程和土水特征曲线进行求解[19]。膨胀土除具有一般黏土的力学行为外，还会因湿度增加而表现出显著的膨胀性，如何将两者相结合并建立膨胀土本构模型是用于后续数值模拟分析的关键。目前，国内外研究者为描述非饱和膨胀土的应力应变关系，对本构模型开展了大量研究，最具代表性之一的是ALONSO 等[20]提出的非饱和膨胀土弹塑性本构模型，该模型可模拟不同应力状态下非饱和膨胀土的反复胀缩变形，能较好地解释膨胀土宏观与微观变形间的相互作用关系，也可用于解决膨胀性非饱和土的变形和渗流耦合的问题[21-22]。但该模型所涉及的参数众多，且微观参数不易获取，完成所有参数测试往往需要很长时间[23]，限制了该模型在工程上的应用。另一个具有代表性的模型就是FREDLUND 等[24]提出的非线性弹性本构模型，该模型可在已知基质吸力的前提下计算相应的应力和位移，也便于数值分析计算，被广泛应用于膨胀土边坡数值模拟计算、膨胀变形预测及膨胀土挡土墙土压力计算[25-26]。此外，还有诸如经验模型[27]、基于功能关系的本构模型[28]等，为膨胀土本构模型的建立提供了参考。综上所述，传统的极限平衡法并不能清晰地描述降雨入渗下的非饱和膨胀土边坡浅层渐进性破坏过程，数值计算方法均将膨胀土视为各向同性材料，膨胀各向异性对膨胀土边坡稳定性的影响机制还有待进一步研究。

考虑膨胀的各向异性，本文作者引入竖向和侧向与基质吸力相关的弹性模量，对非线性弹性本构模型进行修正，并基于室内膨胀试验结果，对模型的可靠性进行验证。基于有限元数值模拟计算软件Geostudio SEEP/W 和ABAQUS，通过二次开发，将非饱和渗流场和本构模型用于后续应力场分析。依托湖北安猇路实体工程，对降雨入渗下非饱和膨胀土边坡进行数值模拟。根据摩尔-库仑强度理论，对比分析膨胀各向异性、膨胀各向同性和不计膨胀这3种工况下的边坡稳定性。

1 本构模型及数值分析方法

1.1 考虑膨胀各向异性的本构模型

为表征非饱和膨胀土在增湿膨胀过程中于竖向和水平向所表现出的膨胀差异，考虑到膨胀土的增湿膨胀是由吸力变化引起，对非饱和土非线性弹性本构模型进行修正。采用竖向和侧向2个与基质吸力相关的弹性模量替代原模型中的对应参数，得到考虑膨胀各向异性的非饱和膨胀土本构模型：

式中：εh和εv分别为水平向和竖向的应变，%；σh和σv分别为水平向和竖向的应力，kPa；ua为孔隙气压力，kPa；uw为孔隙水压力，kPa；(σh-ua)和(σv-ua)分别为水平向和竖向的净法向应力，kPa；μ为泊松比；s为基质吸力，s=ua-uw，kPa；E为与净法向应力相关的弹性模量，MPa；Hv和Hh分别为竖向和水平向与基质吸力相关的弹性模量，MPa。

修正后的模型包含4个弹性参数，即：与净法向应力相关的弹性模量E；与基质吸力相关的弹性模量Hv和Hh；泊松比μ。其中，泊松比μ可在饱和条件下通过K0固结试验[29]测得侧压力系数K0后，由μ=K0/(1+K0)反算得到，见式(2)：

式中：ρd0为试件初始干密度，g·cm-3；λ和N为拟合参数。

3个弹性模量均可通过恒定净法向应力下的膨胀试验测得。假设试件先在无净法向应力条件下增湿膨胀，此时，膨胀应变为εmv；然后，在此基础上施加一定的净法向应力，对应的膨胀应变为εv。在该过程中，吸力始终保持恒定，竖向应变的变化均由净法向应力的增量产生，假设膨胀方向为正，由式(1)可得：

即

式中：εmv为无净法向应力条件下增湿膨胀所对应的竖向膨胀应变，%。

弹性模量是材料本身的属性参数，虽受外界因素(应力状态和边界条件等)影响会产生变化，但始终有一个初始值。但在无净法向应力时，式(4)将因分子分母均为0 而变得无意义。前期研究发现，当基质吸力恒定时，(εmv-εv)随(σv-ua)呈非线性变化，两者的关系可用类似邓肯-张模型的形式进行拟合：

将式(5)代入式(4)可得

式中：a和b为拟合参数，由(εmv-εv)随(σv-ua)的变化曲线拟合得到。可见经此处理后，式(6)在无净法向应力时仍有意义。

而在吸力控制下的膨胀试验时，净法向应力始终保持恒定，由式(1)可得：

式中：εv,VST为恒定净法向应力条件下由可控吸力的竖向膨胀试验(VST)测得的竖向应变，%；εv,LST为恒定净法向应力条件下由可控吸力的侧向膨胀试验(LST)测得的竖向应变，%。

对于初始湿密状态确定的膨胀土，竖向膨胀应变均受到净法向应力和基质吸力的影响而发生变化，故4个弹性参数均随净法向应力和基质吸力的不同而变化。图1所示为由VST测得的侧向力测试值和由修正后的本构模型计算得到的侧向力预测值之间的比较结果，可以发现侧向力预测值与实测值相关性较好，各级净法向应力下侧向力预测值和实测值的相关性系数均大于0.9，说明修正后的本构模型可用于后续非饱和膨胀土边坡的数值模拟分析。

图1 侧向压力测试结果与预测结果对比Fig.1 Comparison of predicted lateral stress and measured lateral stress

1.2 基于ABAQUS平台的UMAT子程序

在ABAQUS 中，需由应变和刚度矩阵计算应力，故式(1)的本构模型可写成如下形式：

式中：{Δσ}为净法向应力增量；{Δε}为应变增量；[D]为净法向应力增量{Δσ}的刚度矩阵；[C]为基质吸力s的柔度矩阵；LT为转置算子矩阵；{I}为单位张量。[D]，[C]和LT的表达式为：

其中：d1=E(1-μ)/[(1+μ)(1-2μ)]；d2=μE/[(1+μ)(1-2μ)]；d3=E/[2(1+μ)]。

考虑到弹性模量E有无净法向应力时的表达式不一致，为避免在参数不连续导致求解不收敛，在有限元分析中采用割线刚度技术进行处理，即在每一次迭代中更新割线刚度直到计算收敛，见式(12)。

式中：Ei为第i次迭代输出的弹性模量E；Ei-1为第i-1次迭代输出的弹性模量E；ξ为衰减系数，本文设为0.9；Efi为由式(6)计算的弹性模量E。

与切线刚度法和Newton-Raphson 法相比，割线刚度法更加简单，且能保证足够高的精度而使迭代收敛。为了提高迭代的收敛性，每次迭代步中的弹性模量E由式(12)计算。另外，相邻2 次迭代步间所产生的累积误差被控制在5%以内，并将其作为判断迭代是否收敛的标准，累积误差Etotali按下式计算：

式中：Eitotal为各节点在第i次迭代的误差。

为了获得计算模型的吸力场，先采用Geostudio SEEP/W 开展饱和-非饱和渗流场模拟，用于计算各个节点的基质吸力。然后，将相应各节点的基质吸力导出并保存为.txt文件。同时，在UMAT 子程序中编写用户自定义场(SDVINI)子程序读取基质吸力的.txt文件，并将其存储在状态变量数组中，以便UMAT 子程序调用。UMAT 子程序计算流程图如图2所示。

图2 UMAT子程序计算流程图Fig.2 Flowchart of UMAT subroutine

2 计算模型与参数

2.1 有限元计算模型

依托工程位于湖北省枝江市安福寺至猇亭区古老背公路工程枝江段(简称安猇路)，沿线约4 km路段处于膨胀土区域。设计边坡坡率为1∶1.5，坡高为8.5 m，在施工过程中，共有9 处膨胀土路堑边坡出现坍塌，左、右幅垮塌边坡总长1.7 km。坡率减小至1∶2.5 后，该边坡仍旧出现垮塌。从现场膨胀土边坡发生垮塌情况看，大多从坡脚处开始滑动，坡中失去土质支撑，承受拉应力，产生较大的裂缝；在坡中发生滑坡一段时间后，坡顶边坡土往下迁移，整个边坡的破坏呈渐进式破坏特征。边坡垮塌后，测量坡顶、坡中、坡脚的裂隙深度和宽度，结果表明坡中的裂隙发育最明显，深度和宽度都比坡顶和坡脚的大，最大裂隙深度为1.5 m，宽度为1.2 cm。对深度为2 m 的膨胀土原状样进行直剪试验，得到有效内摩擦角为17.3°，有效黏聚力为7.8 kPa。据此设膨胀土边坡有限元模型坡高为8.5 m，边坡坡率为1∶1.5和1:2.5。以边坡坡率1∶1.5的边坡为例，其饱和-非饱和渗流和应力场计算模型如图3所示。整个边坡划分为大气干湿循环显著区(区域I)和大气干湿循环非显著区(区域Ⅱ)共2个区域，其中区域I的厚度为3.0 m。监测点位分别布置在坡脚、坡中和坡顶，根据现场量测的裂隙深度，3 点距坡表的竖向距离均设为1.5 m。

图3 膨胀土边坡有限元模型Fig.3 Finite element model of expansive soil slope

采用四边形和三角形混合单元对整个计算几何模型进行网格划分。对于坡率为1∶1.5 的边坡模型，划分后共计819个节点和774个单元；对于坡率为1∶2.5的边坡模型，划分后共计1 210个节点和1 146个单元。膨胀土模型的初始含水率设为18%，对应初始吸力为500 kPa。降雨强度参照文献[17]设置为4×10-7m/s，持续7 d，模拟计算7 d 内非饱和膨胀土边坡渗流场，得到吸力分布随时间的变化。设每小时为1 个时步。对应力场进行模拟时，载荷设置重力方向沿Y方向；模型左、右侧边界受X方向约束，底部受X和Y方向约束。

2.2 计算参数

渗流场计算所需参数包括土水特征曲线和渗透系数，如表1所示。其中，区域I、区域Ⅱ的土水特征曲线分别在0 kPa和50 kPa上覆荷载作用下利用应力相关土水特征曲线测试仪实测[30]，由Fredlund&Xing 模型拟合。区域Ⅱ的饱和渗透系数通过常水头饱和渗透试验测得，区域I的饱和渗透系数则参照文献[17]取值。最后，各区域非饱和渗透系数均由Van Genuchten模型预测得到。

表1 有限元模型水力性质参数Table 1 Hydraulic parameters of FE model

为分析膨胀对非饱和膨胀土边坡稳定性的影响，在应力场分析时，考虑以下3种工况。

1)不计膨胀。应力场计算时，不计吸力减小对弹性模量的软化效应。采用线弹性模型E设为6 MPa，泊松比μ为0.31。

2)膨胀各向同性。弹性模量E按式(6)计算，泊松比μ按式(2)计算。与基质吸力相关的弹性模量则视作各向同性为Hiso，Hiso=(1+μ)/(1-μ)·ds/dεv。

3)膨胀各向异性。弹性模量E按式(6)计算，泊松比μ按式(2)计算。与基质吸力相关的弹性模量视作各向异性，按式(7)计算。泊松比随净法向应力和基质吸力的变化见图4。

在先期研究中，通过饱水条件下的K0固结试验[29]，获得式(2)中的拟合参数λ和N分别为31.790和-9.097。据此，由竖向应变随净法向应力和基质吸力的变化，可得泊松比随净法向应力和基质吸力的变化，如图4所示。

图4 泊松比随净法向应力和基质吸力的变化Fig.4 Variation of Poisson’s ratio with net normal stress and matric suction

为获取弹性模量参数，开展2类吸力控制下恒定净法向应力的膨胀试验，即VST 和LST。VST用于测量E，Hiso和Hv，LST 用于测量Hh。在制备试件时，先通过自制的试件成型盒分层静压，再脱模得到边长为70 mm 的立方体预试件。VST 试件直接采用环刀沿竖向(静压方向)取样得到；LST试件则先将立方体预试件旋转90°，再采用环刀沿预试件侧向(垂直静压方向)取样得到。制备的环刀样直径为61.8 mm，高度为30 mm，经脱模后转移至可开合环刀中并放入非饱和固结仪以开展后续试验，具体试验过程与测试结果见文献[2]。图5所示为考虑膨胀作用时应力场计算所需的弹性模量参数。

图5 弹性模量随净法向应力和基质吸力的变化Fig.5 Variation of elastic modulus with net normal stress and matric suction

3 计算结果及分析

3.1 基质吸力分布分析

通过对非饱和膨胀土边坡历时7 d的降雨进行模拟，得到模型的基质吸力分布及其演变规律。图6所示为区域I内典型的基质吸力沿深度的分布随降雨时长的变化规律。

从图6可见：在降雨条件下，边坡浅、表层土体基质吸力迅速减小；随着降雨历时增加，上部体逐步向深层扩散，导致土体的湿化深度逐渐增加。不同坡率的边坡坡脚、坡中和坡顶监测点的基质吸力随时间的变化如图7所示。

图6 基质吸力沿深度的分布随时间的变化Fig.6 Variation of matric suction distribution along depth with elapsed time

层位的水因自重作用而逐步向下迁移，由表层土

从图7可以发现：在第1天内，因上部雨水还未入渗至1.5 m 深度处，3 处监测点的基质吸力基本保持不变；1 d后，1.5 m深度处的土体开始受到雨水的湿化作用，其基质吸力开始显著减小；4 d后，基质吸力的变化开始逐步趋于平稳，此时，3处监测点的基质吸力均达到最小，但仍未达到饱和状态；降雨7 d后，基质吸力由小到大排列位置依次为坡脚、坡中和坡顶。由于靠近坡顶处径流流速不快，径流量不大，雨水渗透深度有限，故其最终的基质吸力最大；与坡顶相比，坡中和坡脚的径流量较大，径流速度较快，导致土体增湿速度较快，且坡中在增湿后雨水因重力作用也将向坡脚聚集，使得坡脚处的基质吸力比坡中的小。比较不同坡率下的基质吸力历时曲线可以发现：当坡率较小时，坡脚处基质吸力出现明显降幅的时间比陡坡的晚，各监测点的曲线差异也比陡坡的小，说明在较小坡率下，边坡的增湿程度更加均匀，减小边坡坡度可以有效减缓雨水在坡脚处的汇集。

图7 监测点基质吸力历时曲线Fig.7 Variation of matric suction with elapsed time at monitoring points

3.2 应力场分析

分析基质吸力分布发现，3处监测点位的基质吸力衰减程度不一，但衰减幅度均较显著。较大的湿度变化势必会使土体表现出显著的膨胀行为。在不同工况下，通过数值模拟得到的监测点侧向应力时程曲线见图8。图中，AS 表示膨胀各向异性；IS表示膨胀各向同性；NS表示不计膨胀；1.5和2.5分别代表坡率为1∶1.5和1∶2.5。

图8 监测点侧向应力历时曲线Fig.8 Variation of lateral stress with elapsed time at monitoring points

从图8可见：不计膨胀作用时，侧向应力基本不随降雨时长的变化而变化；当考虑膨胀作用时，因基质吸力减小，侧向应力将随降雨历时增加而呈非线性增大，相当于在原侧向应力的基础上增加了1个附加膨胀压力。通过比较各向同性膨胀和各向异性膨胀工况，发现各向同性膨胀下所得侧向应力最大值比各向异性膨胀下的小，将膨胀视为各向同性或不计膨胀将低估膨胀对应力的影响。计算结果表明，相较于膨胀各向异性，膨胀各向同性和不考虑膨胀算得的侧向应力最大值将分别减小20.8%～38.3%和73.9%～88.3%；另外，对于边坡的不同位置，其侧向应力的幅值也不相同，坡脚处最大，坡中次之，坡顶最小，这与基质吸力分布的分析结果相一致。

3.3 膨胀各向异性对边坡稳定性的影响

若不计降雨入渗引起的土容重改变，对于处于一定深度的土体来说，其上覆土自重应力始终恒定。假定土容重为20 kN/m2，可得到3个监测点位的竖向应力均为30 kPa。图9所示为3 种工况下坡脚、坡中和坡顶的应力比K0(即侧压力与竖向应力之比)随降雨历时的变化。

从图9可知：若不计膨胀作用，则应力比始终小于1，说明侧向应力始终小于竖向应力，大主应力为上覆土自重应力；当计入膨胀作用的影响时，膨胀土的基质吸力随着降雨历时增加而减小，膨胀潜势得到释放，受周围土体的约束作用而产生膨胀压力，进而促使侧向压力逐渐增加直至大于上覆土自重应力，导致应力主轴由竖向旋转至侧向；在膨胀各向异性工况下，对于坡率为1∶1.5 的边坡，其坡脚、坡中和坡顶的应力比分别在降雨2，5和6 d时均增大至大于1；而对于坡率1∶2.5的边坡，坡脚、坡中和坡顶则分别在降雨3，4和5 d时才出现应力比均大于1的情况；在膨胀各向同性工况下，坡率为1∶1.5 的边坡坡脚、坡中以及坡率为1∶2.5 的边坡坡脚、坡中和坡顶亦出现了应力比均大于1的情况，只是出现的时间存在差异，这说明坡率对于边坡侧向压力有影响，膨胀作用将使得边坡内应力重分布，并伴随主应力轴旋转。另外，坡脚将先于坡中和坡顶达到一定的应力，上下土层的湿度差异将导致膨胀应力沿深度出现不均匀分布，这使得上下土层间可能出现较大的剪应力，继而造成边坡失稳现象。

图9 监测点应力比K0历时曲线Fig.9 Variation of stress ratio K0 with elapsed time at monitoring points

由前述分析可知：在降雨入渗作用下，非饱和膨胀土边坡土体基质吸力逐渐减小，吸水后有膨胀的趋势，但在侧向约束条件下，这种膨胀势以膨胀力的形式表现出来，导致水平土压力显著增加。在此过程中，非饱和膨胀土的莫尔应力圆的变化如图10所示。

图10 膨胀土增湿过程中的应力状态变化Fig.10 Stress regime change of expansive soil during wetting

对于膨胀土边坡中的某个土单元来说，其上覆土自重应力(σv)保持恒定，而侧向力(σL)则因湿度增加所产生的侧向膨胀力持续增加。当基质吸力仍相对较大时(s=s1)，侧向力是小主应力，此时，莫尔应力圆仍处于抗剪强度破坏包络线之下；当基质吸力减小时(s=s2)，非饱和膨胀土发生增湿膨胀，侧向力增大，莫尔应力圆减小；然而，随着基质吸力的继续降低(s=s3)，侧向膨胀压力继续增大，导致侧向力也在增大，直至大于上覆土自重应力，此时，侧向力成为大主应力，主应力轴从竖直方向旋转到水平方向。在此期间，莫尔应力圆直径不断增加，达到甚至超出抗剪强度破坏包络线。

刘特洪[31]给出了膨胀土沿裂隙面的强度参数推荐值(黏聚力c=5～15 kPa，内摩擦角φ=17°)，由朗肯临界极限土压力理论计算得到被动土压力系数Kp为2.0～3.1。詹良通等[7]通过现场监测也发现雨后测得的应力比与该范围内的数值十分接近。本文选择Kp=2.0 作为非饱和膨胀土边坡发生被动破坏的判定标准。结合图9可以发现：将膨胀视为各向同性及不计膨胀作用时，无论边坡坡率如何，应力比均小于2.0；仅将膨胀考虑为各向异性时，应力比才会出现大于2.0的情况；在降雨入渗4 d和5 d后，坡率为1∶1.5的边坡和坡率为1∶2.5的边坡均于坡脚发生局部被动破坏，继而产生破裂面；坡脚上部的土体失去坡脚土体支撑，在自重作用下承受拉应力，使得破裂面从坡脚开始逐步向上延拓，继而出现渐进式滑坡。因此，膨胀作用导致膨胀土侧向应力显著增加，这是触发边坡失稳的重要原因。

4 结论

1)以非饱和土非线性弹性模型为基础，考虑膨胀的各向异性，引入竖向和侧向与基质吸力相关的弹性模量，提出可以反映膨胀各向异性的非饱和膨胀土本构模型。通过ABAQUS 中的UMAT子程序功能实现本构模型的嵌入，并联合Geostudio SEEP/W的饱和-非饱和渗流模拟，提出了反映膨胀各向异性的湿度场-应力场耦合数值计算方法。

2)在降雨条件下，边坡浅、表层土体基质吸力迅速减小；随着降雨历时增加，上部层位的水因自重作用而逐步向下迁移，由表层土体逐步向深层扩散，导致土体的湿化深度逐渐增加。降雨7 d 后，基质吸力由小到大排列位置依次为坡脚、坡中和坡顶。当坡率较小时，坡面增湿更加均匀。

3)考虑膨胀作用时，因基质吸力减小，侧向应力将随降雨历时增加而呈非线性增大。将膨胀视为各向同性或不计膨胀将低估膨胀对应力水平的影响。计算结果表明，与膨胀各向异性相比，膨胀各向同性和不考虑膨胀算得的侧向应力最大值分别减小20.8%～38.3%和73.9%～88.3%。另外，对于边坡的不同位置，其侧向应力的幅值也不相同，坡脚处最大，坡中次之，坡顶最小。

4)降雨入渗造成膨胀土侧向应力显著增加，并伴随发生主应力轴旋转。坡脚将先于坡中和坡顶达到一定的应力。上下土层的温度差异使得膨胀压力沿深度分布不均，可能导致上下土层间出现较大的剪应力，造成边坡失稳。仅在膨胀各向异性工况下，坡率为1∶1.5的边坡和坡率为1∶2.5的边坡才先后于坡脚发生局部被动破坏。上部土体失去坡脚土体的支撑，在自重作用下承受拉应力，使得破裂面从坡脚开始逐步向上延拓，继而出现渐进式破坏。