引言

边坡动力稳定性数值模拟中，有限元强度折减法逐渐成为解决边坡稳定性的热门方法（黄润秋，2007；涂义亮等，2018），此方法的关键在于选取一个合理的土坡失稳判据，由于土坡破坏机制复杂，使得土坡在极限破坏状态时的失稳判据尚未统一。

目前，有三类判据较为常用，分别是塑性区贯通（等效塑性应变）、数值模拟计算不收敛以及特征点的位移突变或位移不收敛（赵尚毅等，2002），通过上述判据来得到土坡的滑动面位置以及边坡强度储备安全系数。现在的多数学者认为，等效塑性应变贯通是边坡破坏的充分非必要条件（赵尚毅等，2005；吕庆等，2008；龙绪健等，2008），且塑性应变幅值的选取标准带有一定的非确定性等人为因素（裴利剑等，2010）。以有限元计算不收敛作为判据的问题是收敛标准受限于选取的软件，收敛误差也可以调整，作为边坡失稳的判据物理意义不明确，受人为因素的影响（郑宏等，2005；王飞阳，潘泓，2016）。特征点的位移突变或位移不收敛，通常特征点的选取都在坡脚、坡面中点和坡肩处（Zienkiewicz et al，1975；宋二祥，1997），但是滑动带的剪出口可能在坡脚上方，从而使得坡脚处的监测点失去意义，且这种情况只能反应边坡局部状态，边坡可能处于局部失稳，并非整体失稳。史卜涛（2016）则是提出了物质点强度折减法，采用坡顶点竖直方向的位移是否突变作为边坡失稳判据，可见从位移角度分析边坡稳定性的合理性。刘新荣等（2016）从能量角度出发，提出基于能量突变的强度折减法失稳判据。张友利等（2020）、崔笑等（2020）、李志平等（2016）、华成亚等（2016）提出塑性应变能判据应用在边坡动力稳定分析中，并用关键点的残余位移突变判据证明其在边坡动力稳定分析中的可行性和计算精度，在文章中。介玉新（2023）提出相对失稳加速度判别边坡稳定性的方法，指出最小安全系数对应的临界滑动面不一定是最先发生滑动的破坏面，失稳加速度临界滑动面可能最先滑动。刘金龙等（2005）在文献建议联合采用特征点处的位移是否突变和塑性区是否贯通作为边坡的失稳判据，动力作用下边坡失稳判据的选取相对更加复杂。

塑性区反映了土体的状态，当地震作用下边坡土体中出现连续贯通的塑性区时，说明边坡可能发生失稳，但是仅根据塑性区贯通判断边坡是否失稳可能过于保守。因为在实际情况中，即使塑性区贯通，边坡也可能通过内部的应力调整而保持稳定。位移云图是一种基于边坡土体位移变化的评价方法，可以直观地反映出地震作用后土坡的残余位移情况和破裂面形态，该方法能够反映出地震作用下边坡的动态响应。因此，通过塑性区贯通和位移云图的联合评价，能够从两个角度全面考虑边坡的应力应变状态和位移变化情况，从而更快捷准确地判断边坡是否处于失稳状态。本文综合上述判据的特点及差异，提出采用以塑性区贯通为前提，结合土坡呈现出的“彩虹状”位移云图，联合作为土坡动力作用下整体失稳的判据。

1.   彩虹状位移云图失稳判据提出

在地震诱发土坡滑坡的实际工程案例中，发现滑动面与圆弧滑动面相近，通常将其近似为圆弧面来考虑（任永强等，2011；邓东平等，2012；邢洋等，2021）。在理论分析中，各种条分法均是根据瑞典圆弧法（整体圆弧法）发展起来的，瑞典圆弧法原始的求解思想是将整个滑动土体作为一个整体，假设土体均质，滑动面为圆弧滑动面，土体呈刚性转动，在滑面上土体处于极限平衡状态，各力对圆心O的力矩平衡。假定土体单元为不变形刚体，最终建立平衡方程和圆弧滑面力矩平衡方程。当滑坡体沿圆弧滑动面发生相对滑动时，总是沿着所处位置与圆心的这一半径圆弧滑动，首先滑体发生滑动的原因是滑体的下滑力矩M_S超过了抗滑力矩M_R，如公式（1）、（2）所列。

式中：W为滑体重力，d为滑坡体重心与圆心的距离，τ_f为抗剪强度，l为滑动面长度，R为圆弧滑动面半径。

在上述理论和实际工程案例的基础上，本文提出的“彩虹状”位移云图失稳判据是基于地震作用下土坡滑动面为圆弧面的假定。滑坡体以O点为圆心产生整体滑移或转动，如图1所示。

图  1  土坡滑动面示意图

Figure  1.  Schematic diagram of soil slope sliding surface

下载: 全尺寸图片 幻灯片

当滑坡体整体转动θ角度时，滑移面上的点C移动到C′的位置，其位移为：

式中：l_cc′为滑动面点的位移，α为在弧度制下所对的弧数，R为圆弧滑动面的半径。

在本研究的模拟中，考虑了不同的地震动特点，地震动可能会影响位移云图的数值大小和范围，但不会改变彩虹状位移云图的基本特征。基于上述理论推理和数值模拟结果，对结果进行分析后得出边坡失稳的特征，即当边坡整体失稳破坏时的位移云图应满足下面三个条件：

1） 滑坡体内的土体单元相对于边坡内部稳定区沿圆弧滑动面发生相对滑动。边坡内部稳定区与滑坡体在数值模拟的永久位移云图中应呈现出明显的圆弧形界限。

2） 滑坡体位移沿着顺坡向呈现出圆弧形的带状分布。分层明显，类似彩虹状，因此将其定义为彩虹状位移云图。

3） 滑坡体各层位移从坡面指向坡体内部有先增大再突变减小的现象，最大值的层带即为边坡滑动面。

例如，图2分别为相邻两个折减系数对应的土坡地震永久位移响应云图。由图可知，此时滑坡体内部土体单元的位移值小于0.01 m，位移值较小，且同一半径圆弧形条带土体单元的位移量差异较大，并未呈现出彩虹状，表明滑坡体没有沿圆弧滑动面发生整体滑动，因此土坡处于局部破坏状态。

图  2  边坡破坏位移云图

（a）局部破坏位移云图；（b）整体破坏彩虹状位移云图

Figure  2.  Cloud map of slope failure displacement

（a） Localized damage displacement map；（b） Overall damage rainbow displacement map

下载: 全尺寸图片 幻灯片

将边坡继续折减到图2b，位移云图呈现出明显的圆弧形边界，满足上述三个条件（1），可认为图2b中土坡发生整体失稳破坏。

2.   算例分析

2.1   模型参数

本文采用FLAC3D软件进行数值模拟，为研究地震作用下土坡失稳判据，建立边坡模型前缘长40 m，高40 m，坡顶后缘长120 m，高60 m，边坡坡高为20 m，以坡角45°的均质土坡模型为例，同时在土坡模型的坡脚、坡面中点以及坡顶处布设监测点记为P1、P2、P3，网格划分如图3所示。

图  3  边坡模型网格划分示意图

Figure  3.  Schematic diagram of grid division for slope model

下载: 全尺寸图片 幻灯片

本文土坡模型材料参数选取四川地区滑坡体的土体类型及相关数据，土体类型为西南地区粘土，土体物理力学参数如表1所示（曲宏略，2015；工程地质手册，2018）。

表  1  边坡土体参数

Table  1.  Slope Soil Parameters

参数	密度 /（kg·m−3）	体积模量 /MPa	剪切模量 /MPa	抗拉强度 /kPa	粘聚力 /kPa	内摩擦角 /°
数值	1 900	70	37	4	32	24

下载: 导出CSV 

| 显示表格

在模拟地震过程中，对于均质土坡这类简单模型，局部阻尼与实际情况接近（言志信等，2011），且局部阻尼不需要估计系统的自振频率，因此更为简单有效，故本文阻尼的形式为局部阻尼。对于岩土体材料，临界阻尼比取值范围一般为2%—5% （张兆鹏，2018）。本文的研究地区及土体类型与文献（曲宏略，2015）相同，且土体参数相近，因此同样采用5%的临界阻尼比来近似表征土体在地震波传播过程中的阻尼作用，通过式（4）计算可得阻尼系数为0.157 （言志信等，2010）。

式中：ξ为局部阻尼系数，D_s为临界阻尼比，$\text{π} $取3.14。

选取Trinidad地震动，时长21.41 s，峰值0.08 g，计算步数为2 141步，步长0.01 s，加速度、速度及位移时程曲线如图4所示。本文先对地震动进行了滤波处理和基线校正，然后从模型底部竖向输入水平向地震动，为了避免边界的地震波反射，模型周围施加自由场边界条件，底部采用静态边界条件，此时需要将加速度、速度时程转换成应力时程后，再施加到模型底部。在模型底部设置监测点，与输入地震动的加速度时程曲线吻合，说明地震波输入正确。

图  4  Trinidad地震动加速度、速度及位移时程曲线

（a） 加速度时程曲线；（b） 速度时程曲线；（c） 位移时程曲线

Figure  4.  Trinidad Ground Motion Acceleration，Velocity，and Displacement Time History Curve

（a） Acceleration time-course curves；（b） Velocity time-course curves；（c） Displacement time-course curves

下载: 全尺寸图片 幻灯片

2.2   动力有限元强度折减法

本文采用FLAC显式差分法求解微分方程，首先将计算区域离散化为若干单元。在荷载作用下，通过时间步长的有限差分形式计算节点速度和相对位移，进而求出应变和应力。通过不断迭代，求解节点的不平衡力，直至达到平衡状态（陈育民，徐鼎平，2009）。

动力力学平衡方程为：

式中：M为结构的质量矩阵；C为结构的阻尼矩阵；K为结构的刚度矩阵；$ \left\{\ddot{u}\right\} $为结构的加速度列阵，$ \left\{\dot{u}\right\} $为结构的速度列阵，$ \left\{u\right\} $为结构的位移列阵，{F（t）}为结构的节点荷载列阵。

当动荷载为输入地震荷载时，基本力学运动方程可写为：

式中，$ {\ddot{u}}_{g} （ t ） $为输入地震加速度。

地震是一种随时间变化的复杂荷载作用，边坡岩土体在地震作用下很快会进入弹塑性状态（张江伟，2016）。边坡的稳定性安全系数，即边坡的安全储备系数（李德生，徐颖，2012）。在动力作用下边坡的研究方法中，郑颖人等（2005）提出的有限元强度折减动力分析法，在不假定滑面的情况下，既能得出边坡安全系数，又能反映岩土体的动力特性。在不断试算的基础上人为判定边坡是否失稳。在数值分析有限元计算的过程中，强度折减法指通过逐渐提高折减系数F，相应地逐渐减小土体的抗剪强度指标，直至边坡处于极限平衡的临界状态。通过失稳判据来判断此时土坡所处的状态，如果土坡失稳破坏，则减小折减系数F，反之，增大折减系数F。

摩尔—库伦屈服准则对边坡内部滑动面形成机理进行了较详细的说明，不论土体内部应力路径如何变化，只要岩土体试样某一斜截面上的法向正应力和切向剪应力满足摩尔—库伦强度破坏包线方程式（9）或满足主应力形式的摩尔—库伦破坏包线方程式（10），则该试样就会沿着此斜面发生剪切破坏（张坤勇等，2017）。

式中：τ为土体的剪切强度；c为黏聚强度；σ为主应力；φ为内摩擦角。

式中：σ₁为大主应力；σ₃为小主应力。

图5为摩尔库伦破坏包线图。

图  5  摩尔库伦破坏包线图

Figure  5.  Mohr Coulomb failure envelope diagram

下载: 全尺寸图片 幻灯片

式中：χ的实际物理意义为莫尔圆与强度包线之间的最短距离。当χ＞0 时，破坏包线与莫尔圆相离，单元未破坏；当χ≤0 时，破坏包线与莫尔圆相切或相交，单元发生破坏。当潜在滑动面上所有土体单元的应力差χ≤0 时，即为边坡的极限平衡状态。

2.3   土坡地震安全系数计算

通过对模型进行多次折减计算，当折减系数F＝1.25时，如图6a，最终时刻的位移云图并未出现彩虹状，而是边坡下部位移较大，上部位移量较小，处于局部失稳状态；当F＝1.26时，如图6b，最终时刻的位移云图呈现出彩虹状，说明此时边坡处于整体失稳状态，则边坡的强度折减安全系数可取为F＝1.25。

图  6  位移云图

（a） F＝1.25时位移云图；（b） F＝1.26时位移云图

Figure  6.  Displacement Cloud Map

（a） Displacement cloud at F＝1.25；（b） Displacement cloud at F＝1.26

下载: 全尺寸图片 幻灯片

3.   对比其他判据结果

为了验证本文所提出方法的准确性，对比计算了其他不同判据下得到的安全系数，并对其差异性进行阐述。其中，数值软件计算不收敛这一判据对于FLAC3D来说可行性较差，在模拟计算中很难出现不收敛现象，因此在文中并未对此判据进行计算、讨论。

3.1   塑性区贯通判据

图7是边坡模型计算后的塑性区云图，土体单元由于受到剪切或者张拉作用发生剪切破坏和张拉破坏，当破坏单元由坡脚贯通至坡顶时认为边坡失稳。从图7可以看到，当F＝1.23时塑性区未贯通，在F＝1.24时塑性区贯通，故强度折减安全系数F＝1.23。

图  7  塑性区云图

（a） F＝1.23时塑性区云图；（b） F＝1.24时塑性区云图

Figure  7.  Cloud Map of Plastic Zone

（a） Cloud map of the plastic zone at F＝1.23；（b） Cloud map of the plastic zone at F＝1.24

下载: 全尺寸图片 幻灯片

3.2   位移突变判据

位移突变判据通常取坡脚、坡面、坡顶的特征点的最终位移与折减系数的关系曲线变化来判别边坡是否失稳，当折减系数的增大到某一特定值时，该点的最终位移突变，则认为边坡发生失稳破坏，对应的折减系数即为边坡的稳定系数。在研究特征点位移量突变这一判据时，多数采用坡顶的水平位移突变作为特征点位移判据，本文取坡脚、坡面中点以及坡顶的监测点，同时监测特征点水平向与竖直向位移。由图8可以看出，初期随着折减系数的增加，P1，P2，P3的水平向位移和竖直向位移都呈增长趋势。当折减系数增大到1.25时，P2，P3水平和竖向位移同时开始明显增大；当折减系数1.26时，位移量同时发生陡增，位移突变在1.26且继续发展。故边坡的强度折减安全系数可判定1.24—1.26之间。

图  8  强度折减系数与特征点位移量关系曲线

（a） 特征点水平x方向位移与强度折减系数关系曲线；（b） 特征点水平z方向位移与强度折减系数关系曲线

Figure  8.  Relation Curve between Strength reduction Coefficient and Displacement of Characteristic Points

（a） Characteristic point horizontal x-direction displacement versus strength reduction factor curve； （b） Characteristic point horizontal z-direction displacement versus strength reduction factor curve

下载: 全尺寸图片 幻灯片

3.3   位移时程曲线判据

图例中1、2、3分别对应P1，P2，P3，地震作用下边坡监测点水平方向的位移曲线，如图9所示。在折减系数F＝1.24时，监测点的水平向位移随着地震动的进行在0点附近规律性变化浮动，并出现一定永久位移，接着继续对土体参数折减，在折减系数F＝1.25下监测点的位移时程曲线可以看出，位移曲线已经发生了一定程度的偏移，在地震过程中不能再回到0点附近，地震结束后产生了一定的永久残余位移，此时边坡应处于极限平衡或临界失稳状态。当折减系数F＝1.26时，上述现象更加明显，监测点位移在原有位置上产生了较大的偏离，故认为边坡的安全系数在1. 24—1.26之间。

图  9  监测点位移曲线

（a） F＝1.24时监测点位移曲线；（b） F＝1.25时监测点位移曲线；（c） F＝1.26时监测点位移曲线

Figure  9.  Displacement Curve of Monitoring Points

（a） Displacement curve of monitoring point at F＝1.24；（b） Displacement curve of monitoring point at F＝1.25； （c） Displacement curve of monitoring point at F＝1.26

下载: 全尺寸图片 幻灯片

4.   讨论与结论

将四类失稳判据得到的稳定系数列于表2中。

表  2  不同判据下的安全系数

Table  2.  Safety Factors under Different Criteria

判据 类别	塑性区 贯通	特征点 位移突变	位移时程 曲线判据	位移云图 彩虹状
安全系数	1. 23	1. 24—1.26	1. 24—1.26	1. 25

下载: 导出CSV 

| 显示表格

通过表2中对比分析不同判据下边坡的安全系数可以发现，采用特征点位移突变、位移时程曲线和位移云图彩虹状判据得出的安全系数差距不大，说明了本文提出的判据与此前的三类判据具有一致性，也验证了其可靠性。郑颖人等提到的塑性区贯通是边坡土体破坏的必要非充分条件，本文所得结果也可以证明这点，以塑性区贯通为判据得出的安全系数相对较小，与其他判据相比误差为1.6%，出现这种差异的原因可能与网格划分的精密程度与计算精度有关。动力作用下的边坡失稳较为复杂，若采用特征点最终位移突变或者位移时程突变作为判据则需要多次折减系数计算，对于位移监测点是否陡增，更多地依赖于研究者主观上的判定，以此作为判据缺乏说服力。

在采用传统的边坡动力失稳识别方法时，往往存在主观性强、受模拟软件影响的缺陷。本文提出了“彩虹状位移云图”的概念，将边坡滑动面与位移云图相关联，实现了一种新颖且直观的边坡失稳识别方法。通过与传统判据的对比，证明了彩虹状位移云图判据求得边坡安全系数的可行性和准确性，且物理意义明确，操作简便。这使得复杂的数值模拟结果进行解析变得更为简单，提升了工程实践中的效率。

通过上述研究，可得到如下结论：

1） 采用塑性区贯通、特征点最终位移突变、位移时程曲线和彩虹状位移云图四种判据得出的安全系数差距不大，其中采用位移突变判据得到的安全系数具有主观差异性，可取范围标值。

2） 彩虹状位移云图的出现可作为地震作用下土坡整体失稳破坏的判据，物理意义明确，操作简便。同时，也可结合塑性区贯通等其他判据来确认安全系数的准确性。

3） 用特征点位移突变判据时，在选择特征点时应在坡面多选几处，坡脚处的监测点可能会因处于剪出口在其上方而导致结果不准确。

本文的研究对象为简单均质土坡，验证了本文提出的彩虹状位移判据的准确可靠性，对于岩质边坡及多级分层的非均质土坡等还有待进一步研究。