基于强度折减法的坝体边坡稳定性分析

Stability Analysis of Dam Slope Based on Strength Reduction Method

摘要: 土石坝作为水利工程与防洪体系的核心构筑物，其边坡稳定性直接关系到下游区域的生命财产安全与社会经济发展。研究采用有限元分析方法–强度折减法，将土体抗剪强度指标进行折减，以确定边坡达到临界失稳状态时的安全系数。同时，结合达西定律与莫尔–库仑屈服函数，准确描述了渗流场与土体应力应变关系。结果表明，当安全系数超过1.85时，土石坝位移显著增加，出现滑移面，且等效塑性应变集中在顶部，表明此时堤坝进入失稳状态。此研究为土石坝设计、风险预警及加固策略提供了可靠理论依据，有助于提升土石坝安全性与稳定性。

Abstract: As the core structure of water conservancy project and flood control system, the slope stability of earth rock dam is directly related to the life and property safety and social and economic development of the downstream area. In this paper, the finite element analysis method strength reduction method is used to reduce the soil shear strength index to determine the safety factor when the slope reaches the critical instability state. At the same time, combined with Darcy’s law and Mohr Coulomb yield function, the relationship between seepage field and soil stress and strain is accurately described. The results show that when the safety factor exceeds 1.85, the displacement of the earth rock dam increases significantly, and the slip surface appears, and the equivalent plastic strain is concentrated on the top, indicating that the dam enters an unstable state at this time. This study provides a reliable theoretical basis for the design, risk warning and reinforcement strategy of earth rock dams, and helps to improve the safety and stability of earth rock dams.

1. 引言

土石坝作为水利工程与防洪体系的核心构筑物，其边坡稳定性直接关系到下游区域的生命财产安全与社会经济发展。长期以来，受水文地质条件、材料参数不确定性及渗流–应力耦合效应的影响，边坡失稳机理复杂，其安全评估始终是岩土工程领域的重点研究方向。在传统边坡稳定性分析中，极限平衡法(Limit Equilibrium Method, LEM)因其理论简明、计算高效而被广泛应用。该方法通过假设滑裂面形态(如圆弧形、折线形)，基于静力平衡方程计算抗滑力与下滑力的比值(安全系数)。然而，随着工程需求向精细化、动态化发展，极限平衡法的局限性逐渐凸显：其一，其假设的滑裂面形态与真实破坏模式可能存在显著偏差；其二，未能考虑土体应力–应变关系及渐进破坏过程；其三，对复杂地质条件(如非均质材料、渗流场耦合)的适应性不足。因此，亟需一种能够反映材料非线性行为、自动捕捉临界滑裂面且适用于多场耦合分析的稳定性评价方法。

近年来，强度折减法(Strength Reduction Method, SRM)作为一种基于数值模拟的稳定性分析方法，因其独特的理论优势受到广泛关注。强度折减法作为评价边坡稳定性的一种常用数值模拟法，最早由Zienkiewicz [1]于1975年提出，其基本原理是将抗剪强度指标同除以相同的折减系数R，通过改变R值来调整岩土体的抗剪强度从而使边坡达到临界失稳状态，将这个折减系数R值定义为边坡的安全系数Fs。1992年Matsui等人[2]采用Zienkiewicz方法分析了多个边坡的稳定性，并讨论了临界强度折减系数和极限平衡法得到的安全系数的关系，将该法正式命名为“强度折减技术”。随后又有许多学者对该法做了相关研究；如Duncan等[3]、Griffiths等[4]详细论述了将强度折减技术和有限单元法相结合分析边坡的稳定性。Dawson等[5]将强度折减法的结果和上限极限分析法的结果进行对比分析，发现强度折减法的结果略大于极限分析法，并考虑了关联性对强度折减法的影响。Han等人[6]则是对比了不同工况下极限平衡法和强度折减法在滑动面上的差别。国内也有众多学者对强度折减法进行了研究，较早的有宋二祥[7]采用有限元法计算土工结构的安全系数，并在计算中讨论了弧长控制法的应用。随后如郑颖人[8]、唐芬[9]、郑宏[10]等人分别从折减原理、精度控制、屈服准则、失稳判据、三维边坡中应用等方面做了深入的研究，弥补了强度折减法的不足，推动了该法在实际工程中的应用，掀起了强度折减法的研究热潮。张海满[11]运用FLAC^3D软件内置的强度折减功能对土质边坡稳定性算例进行求解分析，结果表明在以上两种水压力条件作用下，使用强度折减法计算所得的边坡安全系数更便捷，且结果可靠。刘航宇[12]采用有限元强度折减法，对一个节理岩质边坡进行了稳定性分析，结果表明，有限元强度折减法可以有效地模拟岩质边坡的破坏过程，考虑了边坡的复杂几何形状、材料非线性、节理裂隙等影响因素，提高了分析的精度和可靠性，为边坡工程设计和施工提供了参考依据。陈诚[13]利用FLAC^3D软件选取两个分区的典型剖面建立边坡的有限元模型，通过对比不同剖面的最大剪切应变增量云图，发现潜在滑动面的位置和范围，从而评估边坡的整体稳定性。通过本研究及时发现边坡的不稳定迹象，采取相应的工程措施，确保边坡的安全稳定。以期为进一步研究露天边坡稳定性提供理论基础和实践指导。陈宁等[14]基于边坡稳定性能的评估问题，介绍有限元强度折减法的基本原理以及在有限元软件上求解边坡安全系数的过程。通过对工程实例计算分析得到最危险的滑动面和安全系数。计算结果与瑞典条分法进行坝坡稳定性分析结果对比表明：两种方法所得的安全系数相差1.47%，滑动面形状、走向和位置一致，并且数值模拟显示的土石坝塑性渐进变化过程符合实际工程塑性变化过程。因此，基于有限元强度折减法对土石坝坝坡等类似工程边坡最终状态稳定性评估分析具有一定的实用性和可靠性。麻荣永等[15]使用ABAQUS有限元软件实现了以强度折减法求解坝坡安全系数，以渗流理论分析了土石坝的渗流稳定，以及求解了考虑流固耦合的土石坝坝坡安全系数。结果表明：有限元强度折减法比极限平衡法更合理、更准确，采用ABAQUS对非饱和土进行稳定渗流分析是合理的、可行的，考虑流固耦合求解土石坝坝坡安全系数更符合工程实际情况，所得结果更加准确。雷艳等[16]对赤金峡水库土石坝开展边坡稳定有限元强度折减分析，对赤金峡水库土石坝河床坝段和岸坡坝段分别开展稳定渗流期、水库水位降落期以及地震作用下的上下游边坡稳定计算，获得坝坡安全系数及塑性剪应变分布等计算结果。基于计算结果评价坝坡的安全稳定，所得结果可以为相似工程提供参考借鉴。

尽管强度折减法在边坡工程中已有较多应用，但其在土石坝稳定性分析中的实践仍面临挑战：一方面，堤坝多由松散填筑材料构成，渗流–应力耦合效应显著，需建立高精度多场耦合模型；另一方面，传统强度折减准则对非饱和土体强度参数的定义尚存争议。针对上述问题，本文以典型土石坝为研究对象，基于COMSOL Multiphysics平台构建渗流–应力耦合数值模型，结合改进的强度折减准则，系统分析边坡的稳定性演化规律。通过对比极限平衡法计算结果及实际工程监测数据，验证强度折减法的可靠性，并揭示渗流场时空分布对安全系数的调控机制，以期为土石坝设计、风险预警及加固策略提供理论依据。

2. 土石坝边坡模拟二维建模

2.1. 有限元强度折减法

所谓强度折减法就是将土体的抗剪强度指标$c\text{'}$ 和$\phi \text{'}$ 用一个折减系数F_s如式(1)和(2)所示的形式进行折减，然后用折减后的抗剪强度指标$c\text{'}$ 和$\phi \text{'}$ 取代原来的抗剪强度指标c和φ在有限元分析安全系数以材料强度的储备方式表达。直到其达到失稳临界状态为止同时得到安全系数F_s。

$c\text{'}=\frac{c}{Fs}$ (1)

$\tan \phi \text{'}=\frac{\tan \phi }{Fs}$ (2)

式中：c——土体的粘聚力；φ——土体的内摩擦角。

有限元强度折减法分析边坡稳定性关键在于根据计算结果判别失稳状态，现有两类失稳判据：一是以力和位移不收敛为标志；二是以广义或等效塑性应变从坡脚到坡顶贯通为标志，第二类是边坡破坏必要条件，第一类合理，本文采用力或位移不收敛临界状态作失稳收敛准则。

2.2. 达西定律及莫尔–库仑屈服函数及相关塑性势

边坡降雨入渗一般需要经历较长的一段时间，边坡土体将经历从干燥到非饱和，再从非饱和状态到饱和状态的转变，其中边坡的物理力学性质会发生转变，一般来说，边坡土体的强度会随着降雨量的增加而逐渐下降，非饱和土体分为固相(颗粒)、液相(水)和气相(空气)，边坡降雨条件下，水的增加将导致土体中气体体积减少，气体的存在将导致水流流速变缓，导致水流运动的路径变长，因此气体的存在将会阻碍降雨入渗的流速，而土体本身的黏滞力也将阻碍水的流动，从而降低土的渗透系数。不论是饱和土还是非饱和土，渗流分析都符合达西定律，达西定律主要定义了饱和土水力梯度与渗流速度的线性关系，即：

$Q=kA\left(h1-h2\right)/L$ (3)

式中：Q为渗流量；A为横截面积；L为渗流长度；$h1-h2$ 为水的高差；k为渗透系数，取决于土的结构和流体的性质与温度。

用v表示单位时间内经过土体横截面积的渗流量，因此达西定律表示为：

$v=\frac{Q}{A}=ki$ (4)

使用平面应变近似在二维模式下对堤坝进行建模。重力和静水压力的影响也包含在内。莫尔–库仑模型的材料属性根据安全系数参数FOS进行参数演化。参数化研究逐渐增加FOS参数，从而降低了每个参数步骤中的土壤强度。对于1.92以上的FOS值，模型不收敛，这表明堤坝失稳。

$F=Q=m\sqrt{J2}+\alpha I1-k$ (5)

其中

$\alpha =\frac{\sin \phi }{3},k=C\cdot \cos \phi$ (6)

参数化内聚力c和内摩擦角为φ：

$C=\frac{c}{Fs},\phi =a\tan \left(\frac{\tan \phi u}{Fs}\right)\left(p<0\right)+a\tan \left(\frac{\tan \phi s}{Fs}\right)\left(p\ge 0\right)$ (7)

式中：c——材料内聚力；φ_u和φ_s分别是非饱和土及饱和土的内摩擦角；p——达西定律给出的孔隙压力。

2.3. 模型建立

本文对某土石坝边坡进行稳定性分析，图1显示某土石坝的横截面。尺寸L1、L2和L3分别为24 m、5 m和24 m，堤坝高度L4为12 m。水位为10 m，可能渗流高度为4 m。堤坝总宽度为L1 + L2 + L3 (53 m)。

Figure 1. Cross section of earth rock dam

图1. 土石坝横截面

该模型中的材料参数以某土石坝边坡为参考取值。其中：杨氏模量E为100 MPa，土壤形变系数即泊松比μ = 0.39，堤坝固体边坡密度ρ_s取2100 kg/m³，堤坝固体边坡孔隙率φ = 0.33，内聚力c = 25 kPa，水密度ρ_w =1 000 kg/m³，安全系数FOS = 1 (安全系数取值范围为1~1.92，安全系数大于1.92时，数值不收敛，可认为堤坝边坡已失稳滑塌)。

(1) 在水库水位至固定土壤域之间设置水体，在水体与堤坝边坡相互作用处(即渗流处)施加扰动模型，即流固耦合发生处。

(2) 流体渗透模型遵从达西定律，流体压力水头和重力作用于堤坝边坡处。将水压作为边界载荷添加到水坝的下游侧。可能的渗流位H_s = 4 m。

(3) 固体边界压力设置在堤坝左侧，下界固定约束，包含重力。土壤塑性屈服准则遵循莫尔–库仑准则。外部应力为孔隙压力和重力。

对模型进行网格划分，利用有限元分析得到完整网格计算模型，网格划分采用自由三角形网格，模型共有节点44,280个，网格包含9686个域单元和276个边界层单元，如图2(b)所示。

(a) 模型几何示意图

(b) 网格划分

Figure 2. Model schematic diagram and grid division diagram

图2. 模型示意图及网格划分图

3. 结果与分析

3.1. 压力水头等值线分析

土石坝压力水头等值线及零压力水头线如图3所示，压力水头以堤坝边坡左侧受力面为受压浸水区，在浸没壁上压力水头从堤坝顶部向下延伸至下部固定域。在浸没壁上，它从0 m到10 m不等，而在渗流面上为0 m。正压力水头表示正孔隙压力，表示饱和土，而非饱和土用零压力水头表示。图中的零压力水头线是将饱和土与非饱和土分开的潜水面位置。

Figure 3. Pressure head contour

图3. 压力水头等值线

3.2. 安全系数与最大位移的关系

Figure 4. Relation curve between maximum displacement and safety factor

图4. 最大位移与安全系数的关系曲线图

图4和图5分别显示了安全系数(FOS)与堤坝最大位移之间的关系，以及土石坝失稳时的位移大小分析结果。图4中，横轴表示安全系数，纵轴表示最大位移(单位为毫米)。在安全系数为1到1.85之间，最大位移几乎为零，表明在这个范围内堤坝是稳定的。然而，当安全系数超过1.85时，最大位移迅速增加，表明堤坝开始失稳。当安全系数为1.85时，最大位移约为5毫米；当安全系数增加到1.9时，最大位移迅速增加到110毫米。最大位移在FOS = 1.9附近显著增加，这表明堤坝开始失稳。土石坝边坡临界失稳时，土石坝右侧下部土体最先开始发生滑移，由于孔隙水压力的作用，下部土颗粒发生位移，导致渗流破坏进一步扩大，土石坝上部土体缺少支撑，开始滑移，堤坝发生失稳滑塌[17]。

Figure 5. Diagram of critical instability displacement of slope

图5. 边坡临界失稳位移大小图

3.3. 等效塑性应变以及滑移面分析

Figure 6. Equivalent plastic strain diagram

图6. 等效塑性应变图

边坡失稳时滑弧及滑移面由图6和图7可知，失稳前的等效塑性应变表明了破坏机理。滑移面箭头表示土颗粒的位移方向，说明了边坡土体滑移现象，由于边坡下边界为固定约束，因此边坡右下角土体不会出现滑移。滑移面起点为左侧边坡渗流处，沿潜水面斜向下滑移。这种等效塑性应变的分布表明边坡右下角的土壤已经进入塑性变形阶段，可能导致滑移。滑移面主要分布在土石坝顶部，表明土石坝顶部是潜在的失稳区域。滑移面在土石坝顶部的深度约为10米，沿着土石坝的坡面延伸。这种滑移面的分布表明土石坝顶部的土壤受到较大的剪切力导致滑移。

Figure 7. Sliding surface of earth rock dam in instability

图7. 土石坝失稳时滑移面

4. 结论

(1) 基于莫尔–库仑强度准则，边坡稳定性分析中材料参数(黏聚力、内摩擦角)与安全系数(FOS)呈负相关演化规律。在库区堤坝稳定性数值模拟中，定义边坡滑移位移为FOS的函数变量时发现：FOS的提升会加剧滑移面位移响应，进而引发土体参数劣化与抗剪强度衰减，最终触发堤坝失稳临界状态。

(2) 采用二维平面应变模型对堤坝进行稳定性预测时，可通过后处理技术实现滑移路径的三维空间重构。数值计算收敛性可作为失稳判据的量化指标，研究数据表明：当FOS > 1.92时，弹塑性本构模型的计算稳定性显著降低，此时对应堤坝结构破坏的起始阈值。

(3) 在特定荷载工况下，堤坝稳定性与FOS存在临界突变效应(本研究中临界值FOS = 1.92)。该现象源于塑性应变累积与剪切强度非线性衰减的耦合作用机制。位移–FOS关系曲线显示：当FOS接近1.9时，特征点位移量呈现数量级跃升，此突变特征可作为失稳预警的量化依据。