线性与非线性强度的土石坝坝坡稳定分析下限法

周建烽，王均星，陈 炜，罗贝尔

（1.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；2.福建省水利规划院，福建 福州 350001）

1 引 言

目前，土石坝稳定分析在工程上广泛采用的是刚体极限平衡法[1－4]。刚体极限平衡法概念清晰、计算简便，工程实践经验比较丰富。但它假定失稳块体为整体移动的刚体，只考虑滑移体上力的平衡，不考虑结构的屈服，不能确定相应的应力和位移分布。为求改进，国内外众多学者都在寻求一种更有效的稳定分析方法，其中最受关注的是有限元法和塑性极限分析方法。在有限元法中，一般采用有限元强度储备系数法，它可以不作任何假定，能够考虑土体的本构关系，能够模拟边坡的破坏过程及其滑移面形状，并且可以很方便地考虑土体和锚固、支护的共同作用。采用有限元强度储备系数法进行边坡稳定分析是近年来的热点，国内外许多学者都进行了有益的探索[5－9]，尤其是可以利用它来进行边坡失稳的仿真计算。

塑性极限分析方法避开弹塑性分析的全过程，直接研究结构的极限状态，求解结构的极限荷载，为评价结构的安全度提供必要的数据，其理论严谨，在土工问题极限承载能力上是一种十分有效的方法[10]。自引入有限元法解决了构造静力场和机动场的难题之后，结合数学规划手段，已成功地用于计算地基承载力及边坡稳定分析等方面。

本文采用塑性极限分析方法进行土石坝的坝坡稳定分析，在Sloan[11]的工作基础上，考虑了渗流作用和地震荷载，推导了基于非线性规划的土石坝有限元塑性极限分析下限法模型，并采用迭代算法对堆石、砂土等无黏聚土进行非线性强度指标的坝坡稳定计算。最后，对几个典型土坡和具体的土石坝工程进行了计算，并与多种方法的分析结果相比较。

2 塑性极限分析的下限原理

在物体及其边界上，满足下列条件的应力场称为静力容许应力场：

（3）在力边界上，荷载与真实荷载相同。

静力容许应力场有多个，每个静力容许应力场对应着一个外荷载。因为结构并未破坏，因而这个荷载一定比极限荷载小。下限定理可以表述为：在所有的与静力容许应力场对应的荷载中，最大的荷载为极限荷载，与最大荷载相应的静力容许应力场为结构的极限状态应力场。

3 土石坝坝坡稳定的有限元极限分析下限法模型

土石坝坝坡稳定的有限元塑性极限分析下限法非线性规划数学模型由三部分组成：①非线性约束条件，即屈服条件；②线性约束条件：包括平衡条件、应力间断条件和边界条件；③目标函数，即安全系数。本文采用线性三角形单元离散结构体，每个单元都需满足上述的约束条件与目标函数。

3.1 结构的离散

为构造静力许可应力场，结构物采用不共节点的三角形线性单元离散，以节点应力为未知量。考虑渗流作用后，每个节点都增加一个水头变量Hi。

图1 下限分析法单元离散图Fig.1 Element discrete of lower bound limit analysis

在下限法中，单元内任意一点的应力分量、水头可表示为节点应力分量、节点水头的线性函数，即

式中：A为三角形单元的面积，ξi、ηi、ζi的具体表达式参见一般有限元著作。

3.2 平衡约束条件

土石坝稳定计算要考虑多种荷载的影响，其中最主要的外荷载为孔隙水压力和地震荷载。

对于孔隙水压力的考虑，采用有效应力原理，根据水土分算原则，即把土骨架当作有限元极限分析法中的研究对象，把孔隙水压力当作是一种外力作用于土骨架上。本文在考虑渗流作用时，将其视做体积力考虑并表示为

式中：γw为水的重度；J为水力坡降。

可以把Fp分解为两个分量：

式中：Fpx、Fpy分别为渗流作用力在x、y 方向的分量； Jx、 Jy分别为水力坡降在x、y 方向的分量。

土石坝用动力法作抗震计算时，由于坝料的非线性特性使得计算仍然十分复杂。迄今为止，对于土石坝的地震稳定性和永久变形的动力方法和准则在工程界仍然没有共同的看法。我国的水工建筑物抗震规范[12]规定土石坝抗震稳定计算可采用拟静力法。拟静力法将地震荷载作为静荷载来检验结构的强度，在静力法的形式中，纳入了动力分析法的内容，它具有简便、实用的优点。因此，本文在有限元极限分析下限法中，采用了拟静力法来考虑地震荷载。

下限法平衡方程表示为

式中：kh、kv分别为水平和垂直地震加速度代表值；γ′为土体重度（浸润线以下取浮重度，浸润线以上取天然重度）；γ为包含孔隙水在内的土体重度，应力拉正、压负。

平衡条件写成矩阵形式如下：

式中：上标e表示单元号；且有

对于其他线性约束条件即应力间断条件和边界条件的建立，可以参考文献[13]，本文不再赘述。

3.3 屈服约束条件

对于一般受力下的岩土材料，通常采用Mohr-Coulomb屈服准则。在有限元极限分析下限法模型中，对于平面应变问题，每个单元节点应满足屈服条件，可表示为

式中：i=1，2，3，e=1，2，…，NE，NE为总单元数，分别为第e 个单元的3个应力分量、材料的黏聚力和内摩擦角。

3.4 目标函数及非线性数学规划模型

本文采用强度储备系数来定义安全系数，即

式中：φ0、c0为进行强度折减以后的强度参数。

式（13）可简写为

将强度储备系数设为目标函数的下限法数学规划模型为

对上述形成的土石坝坝坡稳定分析下限法的非线性数学规划模型求解，本文采用专业优化软件Lingo8.0，通过非线性规划求解最终可以获得土石坝坝坡抗滑安全系数值及相应的静力许可应力场。

4 非线性强度指标的坝坡稳定计算

高土石坝一般多采用堆石料，对此类材料在高土石坝的高应力水平条件下采用线性强度指标已不能满足要求，有必要采用非线性的强度指标。

本文在考虑材料强度的非线性时，采用Duncan对数模式为

式中：φ为土体滑动摩擦角（°）；φ0为一个大气压下的摩擦角； Δφ为增加一个对数周期下φ 的减小值；σ1为土体的大主应力值（本文以拉应力为正）；Pa为大气压力。

该模式在我国有十分丰富的应用经验，碾压式土石坝设计规范规定对于粗粒料非线性抗剪强度指标采用该模式。

在下限法模型中，本文采用了迭代算法进行非线性强度指标的坝坡稳定计算，具体步骤如下：

（1）拟定一个初始内摩擦角 φ0，一般取φ0。

（2）用本文第3节中的方法进行计算，得到第j 次迭代的安全系数 kj及应力场{σj}，并获得第一主应力

（3）由式（16）计算出各单元内摩擦角φe。

（4）利用步骤（3）的结果重复步骤（2），得到第j+1次迭代的安全系数kj+1及应力场{σj+1}，验算前、后两次安全系数是否满足，ε为收敛精度。

（5）如果满足收敛精度，则停止迭代，得到的安全系数kj+1即为所求的安全系数。如不满足则重复步骤（2）～（4）。

5 算例分析

5.1 含软弱夹层的算例

图2所示为一含有软弱夹层的土坡，Fredlund[13]、Kim[14]等分别采用刚体极限平衡法及基于线性规划的极限分析法对此算例进行过分析。该土坡顶宽L=18.3 m，高H=12.2 m，水头高度HW=6.1 m，坡比为1:2，均质土体材料参数为c′=28.76 kPa，γ=18.89 kN/m3，软弱夹层材料参数c′=0 kPa，φ′=10°。

本文采用3种不同密度的网格（74个单元、156个单元、238个单元，如图3所示）进行计算，计算结果见表1。由本文方法还可以得到极限状态下的应力场。由于篇幅原因，此处只给出密网格工况3、4的第2主应力等值线图，详见图4、5。由主应力等值线图可看出，在孔隙水压力和土体自重作用下，主应力大小随着深度加深而线性增加，符合一般自重应力场的分布规律。

从表1的计算结果可知，安全系数随着网格的加密而增加，并且当单元数达一定程度时，增加网格数量提高的精度并不明显。本文方法在网格密到一定程度时计算的结果比瑞典法高2%～5%左右，比简化毕肖普法和摩根斯顿-普赖斯法小3%～5%左右。从塑性极限分析的观点来看，刚体极限平衡法考虑了力的平衡条件，在滑动面上满足屈服条件，其解答可归为下限法的范畴，但并不能判断其他部位是否满足屈服条件，而且瑞典法就连平衡条件也没有完全满足，因而不是严格的下限解；另外，它假设了滑动面，即假设了结构破坏的一种机构，并且从一系列破坏机构中寻找安全系数的最小值作为最终解，其解答也可归为上限法的范畴，但它并没有考虑变形协调关系，事实上也不是一种完全的上限解，因而由极限平衡法得到的结果既非下限解也非上限解。

由此可以说明，本文方法得到的结果为严格的下限解，是偏于安全的。本文结果比Kim等[15]计算的下限解略高，这是因为Kim在对屈服函数的线性化时采用内接正多边形对屈服圆进行逼近，根据下限定理得到的解必然比本文方法小。

图2 含有软弱夹层的土坡示意图Fig.2 Sketch of slope with a weak interlayer

图3 有限元网格划分Fig.3 Meshes of finite elements

表1 多种方法的安全系数Table 1 Factor of safety for slopes obtained with different methods

图4 工况3第2主应力图（单位:kPa）Fig.4 Contours of second principal stress in slope at condition 3 (unit:kPa)

图5 工况4第2主应力图（单位：kPa）Fig.5 Contours of second principal stress in slope at condition 4 (unit:kPa)

5.2 线性和非线性强度指标的土坝坝坡稳定计算

某土坝坝坡坡角为α，坡比为1:2，坝料的线性强度参数c=0，φ=45°，非线性强度参数为c=0，φ0=55°，Δφ=10°。本文计算不同 γ H值下的安全系数，以获得安全系数随 γ H变化规律，计算结果见表2。图6、7分别为 γ H=4 000 kPa情况下的第1和第2主应力图（本文以拉应力为正）。由图可看出，主应力（有效应力）等值线分布良好，符合一般规律。

由表2中结果可知：

（1）当采用线性强度指标计算时，本文方法与瑞典圆弧法和简化毕肖普法计算出的安全系数极为接近，皆近似等于tan φ/tanα，由此可以看出，本文方法计算结果是可靠的。

（2）采用非线性强度指标的计算结果两者相差较大，达15%左右，并且在 γ H较大时，本文算出的结果比后者要小；在 γ H较小时，本文算出的结果反而大。究其原因，刚体极限平衡法无论是进行线性还是非线性计算，仅能计算滑动面上的应力，只满足滑动面上的屈服条件。其计算结果既不是下限解也不是上限解。

表2 线性和非线性强度指标下的坝坡安全系数Table 2 Dam slope factors of safety with linear and nonlinear strength indexes

图6 第一主应力图（单位：kPa）Fig.6 Contours of the first principal stress in slope (unit:kPa)

图7 第二主应力图（单位：kPa）Fig.7 Contours of the second principal stress in slope (unit:kPa)

5.3 土石坝稳定分析算例

土石坝规范[16]规定土石坝稳定计算应分别对施工期（包括竣工期）、稳定渗流期、水库水位降落期和正常运用遇地震4种工况进行计算。由于本文研究目的在于验证塑性极限分析这一方法在土石坝稳定分析中的可行性，所以只分别对心墙坝的稳定渗流期遇地震工况、斜心墙坝的正常蓄水位遇地震工况进行计算。至于其他工况，不同之处只是在坡外水水头、材料参数和孔隙水压力取值上的差别。

5.3.1 心墙坝算例

某心墙坝坝体断面图和材料分区如图8所示，坝体力学参数见表3，坝高为32.5 m，上游水位为27.5 m，下游无水，水的重度取为9.81 kN/m3，地震烈度为Ⅷ度。

表3 土体力学参数Table 3 Mechanical parameters of soils

计算模型中共划分了738个单元，共有2 214个不共用节点，1 069条公共边（见图9）。由下限法计算得到的坝坡抗滑安全系数见表4（本文方法考虑的是土石坝上下游整体，不是某个方向的滑动），单元主应力矢量见图10。

图8 心墙坝断面示意图Fig.8 Sketch of core wall dam section

图9 心墙坝有限元网格划分Fig.9 Meshes of finite element of core wall dam

表4 心墙坝坝坡抗滑安全系数Table 4 Factors of safety for core wall dam

图10 心墙坝主应力矢量图Fig.10 Vectogram of principal stresses of core wall dam

5.3.2 斜心墙坝算例

某斜心墙坝坝体断面图和材料分区如图11所示，坝体及坝基材料物理力学性质指标见表5。水的重度取为9.81 kN/m3，地震烈度为Ⅵ度。

模型共划分了901个单元，共有2 703个不共用节点、1 311条公共边（见图12）。由下限法计算得到的坝坡抗滑安全系数见表6，单元主应力矢量见图13。

图11 斜心墙坝断面示意图Fig.11 Sketch of inclined core wall dam section

表5 坝体及坝基物理力学性质指标表Table 5 Physico-mechanical parameters of dam body and foundation

图12 斜心墙坝有限元网格示意图Fig.12 Meshes of finite element of inclined core wall dam

表6 斜心墙坝坝坡抗滑安全系数Table 6 Factors of safety for inclined core wall dam

图13 斜心墙坝主应力矢量图Fig.13 Vectogram of principal stresses of inclined core wall dam

由心墙和斜心墙土石坝计算结果可知，由下限法算出的坝坡抗滑安全系数结果与瑞典圆弧法和简化毕肖普法基本相符，主应力分布规律与常规有限元法的结果基本一致，自重应力为主要分布，符合一般的应力分布规律。

6 结 论

（1）本文在Sloan的工作基础上，将有限元塑性极限分析下限法用于土石坝坝坡稳定分析，理论严格且物理意义明确，可求得坝坡抗滑安全系数的下限解及相应的静力许可应力场。通过几个典型土坡和具体的土石坝工程的分析发现，本文方法的结果与极限平衡法相符，并且无论是使用线性强度指标还是非线性强度指标进行计算，都能获得较好的结果，证明了本文方法的正确性。

（2）刚体极限平衡法在一定意义上满足了静力平衡条件，假定了破坏机制，一定程度上满足了上下限定理的一些条件，但并不满足完全屈服条件和变形协调条件，其计算结果既不是下限解也不是上限解。

（3）本文方法具有以下优点：①不用假定滑动面，避免了由于滑动面位置和形状带来的计算误差，并在解决了规划问题的前提下，可以对三维问题进行分析；②能够较方便地处理复杂几何形状和非均质材料问题；③回避了结构材料复杂的本构关系；④可以考虑各种外荷载。

（4）目前，塑性极限分析在水工领域应用还比较少。为发挥其自身优势使其今后能够在水工方面得到推广，今后还需做的工作是：①通过对大量的实际工程进行计算分析，与现行的各种计算方法进行比较后，确定与规范相对应的下限法坝坡抗滑安全系数的取值范围；②由于求解是结合优化方法进行的，寻求一种有效的求解大规模稀疏矩阵的算法和求解器是必要的；③结合上限法构造一套力学概念清晰、计算方法先进、计算结果有效全面的土坝坝坡稳定分析方法。

[1]BISHOP A W.The use of the slip circle in stability analysis of slope[J].Geotechnique，1955，5(1):7－17.

[2]MORGENSTERN N R，PRICE V E.The analysis of the stability of general slip surfaces[J].Geotechnique，1965，15(1):79－93.

[3]SPENCER E.A method of analysis of the stability of embankments assuming parallel inter-slice forces[J].Geotechnique，1967，17(1):11－26.

[4]SARMA S K.Stability analysis of embankments and slopes[J].Geotechnique，1973，23(3):423－433.

[5]TAMOTSU，MASTUI，SAN KA CHING.Finite element slope stability analysis by shear strength reduction technique[J].Soils and Foundations，1992，32(1):59－70.

[6]连镇营，韩国城，孔宪京.强度折减有限元法研究开挖边坡的稳定性[J].岩土工程学报，2001，23(4):407－411.LIAN Zhen-ying，HAN Guo-cheng，KONG Xian-jing.Stability analysis of excavation by strength reduction FEM [J].Chinese Journal of Geotechnical Engineering，2001，23(4):407－411.

[7]郑颖人，赵尚毅.有限元强度折减法在土坡与岩坡中的应用[J].岩石力学与工程学报，2004，23(19):3381－3388.ZHENG Ying-ren，ZHAO Shang-yi.Application of strength reduction FEM to soil and rock slope[J].Chinese Journal of Rock Mechanics and Engineering，2004，23(19):3381－3388.

[8]陈国庆，黄润秋，周辉，等.边坡渐进破坏的动态强度折减法研究[J].岩土力学，2013，34(4):1140－1146.CHEN Guo-qing，HUANG Run-qiu，ZHOU Hui，et al.Research on progressive failure for slope using dynamic strength reduction method[J].Rock and Soil Mechanics，2013，34(4):1140－1146.

[9]李小春，袁维，白冰，等.基于局部强度折减法的边坡多滑面分析方法及应用研究[J].岩土力学，2014，35(3):847－854.LI Xiao-chun，YUAN Wei，BAI Bing，et al.Analytic approach of slope multi-slip surfaces based on local strength reduction method and its application[J].Rock and Soil Mechanics，2014，35(1):847－854.

[10]陈惠发，詹世斌.极限分析与土体塑性[M].北京:人民交通出版社，1995.

[11]SLOAN S W.Lower bound limit analysis using finite elements and linear programming[J].International Journal for Numerical and Analytical Methods in Geomechanics，1988，12(1):61－77.

[12]中华人民共和国行业标准编写组.SL203－97水工建筑物抗震设计规范[S].北京:中国水利水电出版社，1997.

[13]王汉辉，王均星，王开治.边坡稳定的有限元塑性极限分析[J].岩土力学，2003，24(5):733－738.WANG Han-hui，WANG Jun-xing，WANG Kai-zhi.Plastic limit analysis of slope stability using finite element[J].Rock and Soil Mechanics，2003，24(5):733－738.

[14]FREDLUND D G，KRAHN J.Comparison of slope stability methods of analysis[J].Canadian Geotechnical Journal，1977，14(3):429－439.

[15]KIM J，SALGADO R，LEE J.Stability analysis of complex soil slopes using limit analysis[J].Journal of Geotechnical and Geoenvironmental Engineering，2002，128(7):546－557.

[16]中华人民共和国行业标准编写组.SL274－2001碾压式土石坝设计规范[S].北京:中国水利水电出版社，2002.