基于ANSYS的有限元强度折减法确定边坡安全系数

基于 ANSYS的有限元强度折减法确定边坡安全系数摘要：本文基于ANSYS，采用D-P外角点外接圆屈服准则对国内某矿区边坡进行稳定性计算分析，通过不断对边坡强度参数黏聚力和内摩擦角进行折减，直到软件计算不收敛为止，其折减的倍数即为边坡稳定安全系数。

计算结果显示，利用ANSYS自带D-P本构模型计算得到的边坡安全系数远大于极限平衡法计算得到的边坡安全系数。

最后应用不同屈服准则安全系数的转换关系得到该边坡平面应变下与M-C匹配的D-P准则的安全系数，并与极限平衡法结果对比，吻合较好。

据此得出结论：在估算边坡安全系数方面，采用有限元强度折减法是一种值得信赖的方法，但计算中采用理想弹塑性材料模型时，屈服准则的选择会对边坡安全系数的计算产生较大影响。

关键词：有限元；强度折减；屈服准则；边坡稳定；安全系数1 引言目前，边坡稳定性分析发方法较多，主要有定性分析法（图解法、类比法）、定量分析法（极限平衡法、数值分析法）、非确定分析法（模糊分析评判法、可靠性分析法）。

而对于边坡安全系数，许多学者大多用定量分析法[1-2]。

传统的极限平衡法首先要确定一个潜在的滑动面，基于一系列简化假定后，由力系平衡或能量守恒求得滑动面的安全系数，用它作为评价边坡安全性的指标。

这些方法有瑞典条分法、简化毕肖普法、简布法、不平衡系数传递法等。

这些方法的基本出发点是一样的，即刚塑性假定，不同之处在于对条间力所作的假定不同。

由于这些假定的物理意义不一样，因此它们所能满足的平衡条件也不相同，计算步骤有繁有简，为了检验所列的各方法和其他边坡稳定性分析方法的精确性，许多学者在过去几十年里从不同角度做了大量研究并进行了系统总结[3-7]。

传统的极限平衡法由于没有考虑土体内部应力与应变的关系，故无法模拟分析土体发生变形甚至破坏的过程。

随着计算机技术的发展，数值计算方法在边坡稳定分析中得到了广泛的应用。

其最大的优点是求解安全系数时，不需要假定滑移面的形状和位置，也无需进行条分，可以分析任何形状的几何体，不但能进行线性分析还可进行非线性分析。

它不仅考虑了土体的平衡条件，而且考虑了材料的应力与应变关系，使得结果更趋于准确。

目前用于边坡稳定性分析的数值方法主要有：有限单元法、边界单元法、有限差分法、离散单元法和块体系统不连续变形分析法等。

但究其实用性、灵活性和应用的广泛性而言，有限单元法是数值模拟方法中用于岩土力学分析的有效技术之一。

2 有限元强度折减法2.1 强度折减法原理所谓强度折减法[8]就是将土体的抗剪强度指标c和φ，用一个折减系数Fs进行折减，然后用折减后的抗剪强度指标cm 和φm取代原来的抗剪强度指标c和φ，不断折减，直到边坡出现破坏为止，同时得到安全系数Fs。

首先设定一个折减系数Fs ，然后根据折减系数原理，得到cm和φm，将折减后的参数输入到ANSYS的D-P本构方程中，进行有限元计算。

如果计算收敛，则坡体处于稳定状态，不断地增加Fs ，至到计算发散为止，此时的折减系数Fs即为边坡的稳定安全系数，此时的滑移面即为实际的滑移面。

2.2 安全系数的定义传统的边坡极限平衡方法采用Mohr-Coulomb屈服准则，安全系数定义为滑动面的抗剪强度与滑面上实际剪力的比值，用公式表示为式(1)，将式(1)两边同时除以Fs，则式(1)变为式(2)。

(1)(2)其中：cm =c/Fs；tanφm=(tanφ)/Fs。

式(2)表明当抗剪强度参数折减Fs后，坡体达到极限平衡状态，可以看出有限元强度折减法在本质上与传统方法是一致的。

2.3 屈服准则和流动法则选取边坡稳定性分析主要关心的是力和强度之间的关系，对于本构关系的选择不必十分严格[9]，本文计算采用理想弹塑性模型。

但屈服准则的选取对于安全系数的计算影响较大，目前岩土工程中广泛采用Mohr-Coulomb屈服准则，引入应力罗德角参数M-C屈服准则可以表示为：(3)式中，I1为应力张量第一不变量；J2为应力偏量第二不变量；为应力罗德角。

该准则在三维应力空间中由六个分段函数构成，在π平面上的屈服曲线为不规则六角形(见图1)，数值计算极为不便[10]。

因此学者陆续提出D-P系列修正模型(见表1)。

表1的模型有一个通用的表达式：(4)此式是1952年由Drucker-Prager提出的，由于该屈服准则在π平面上的屈服曲线为圆形(见图1），目前国际上流行的大型有限元软件大都采用了Drucker-Prager准则。

可以看出当式(3)中的为不同常数时，屈服函数简化为特定条件下的D-P准则[11]。

D-P系列准则是M-C屈服准则的特殊情况，有关参数参见表1。

D-P系列准则能否正确使用取决与岩土体不同的应力状态，在单向压缩及常规三轴压缩等满足σ1= σ2>σ3条件的应力状态下，DP1 是与 M-C 准则匹配的；在单向拉伸及常规三轴拉伸等满足σ1> σ2=σ3条件的应力状态下，DP4 与 M-C 准则匹配；在平面应变条件的关联流动法则下(膨胀角ψ =φ)，DP5 与 M-C 准则匹配；而在非关联流动法则下(膨胀角ψ =0)，DP3 与 M-C 准则匹配；在实际分析中与 M-C 准则匹配的 DP2 所对应的应力状态没有明显的物理意义，它只是一种与 M-C 准则近似等效的准则，与 M-C 准则的差别程度取决于实际的材料参数与应力状态，计算精度随着实际的材料参数与应力状态而变化。

图1 M-C准则与D-P类准则在π平面上的曲线Fig.1 The curves of M-C and D-P criteria on the π plane 目前美国大型有限元软件ANSYS采用的是Mohr-Coulomb不等角六边形外接圆D-P屈服准则(D-P1)，研究表明评价边坡稳定性等问题该准则偏于不安全[11]，为了和传统工程中Mohr-Coulomb准则下安全系数接轨，郑颖人[9]提出对于边坡的稳定性分析，平面应变状态下的强度问题可采用用D-P3或D-P5，三维空问题下可采用在π平面上与Mohr-Coulomb不等角六边形面积相等的DP2准则。

表1 常用屈服准则Table 1 Commonly used yield criteria材料参数屈服面位置号(相对于M-C屈服准则)αк-CM-C 屈服面P1外角点外接D-P 圆P2M-C 等面积D-P 圆P3平面应变M-C 匹配圆P4内角点外接D-P 圆P5内切D-P 圆本文基于ANSYS 分析结果，得到采用DP1屈服准的边坡安全系数，利用同济大学钟才根[12]研究成果，通过安全系数转换关系式进行转换，得到DP3屈服准则下的边坡安全系数，不需要进行二次开发。

选择由于在 ANSYS 分析中，对于 D-P 屈服准则，需要输入材料的c 、φ 和ψ 参数，其中ψ 为材料的剪胀角。

根据本文选择的屈服准则，自动取非关联流动法则下的剪胀角ψ=0。

2.4 边坡失稳判据采用有限元分析采用有限元强度折减法分析边坡稳定性的一个关键问题是如何根据有限元计算结果来判别边坡是否处于整体破坏状态。

目前,土体破坏的标准有如下几种[9]：①以有限元静力平衡计算不收敛作为边坡整体失稳的标志；②以塑性区(或者等效塑性应变)从坡脚到坡顶贯通作为边坡整体失稳的标志；③土体破坏标志应当是滑动土体无限移动，此时土体滑移面上应变和位移发生突变且无限发展。

根据郑颖人等人的研究，上述土体破坏3种标准有如下关系：土体滑动面塑性区贯通是土体破坏的必要条件，但不是充分条件。

土体整体破坏的标志应是滑体出现无限移动，此时滑移面上的应变或者位移出现突变，因此，这种突变可作为破坏的标志。

此外有限元计算会同时出现计算不收敛。

可见，上述①③两种判据是一致的。

因而可将有限元数值计算是否收敛或者滑面上节点塑性应变和位移突变作为土体破坏的依据。

3 实例分析3.1 实例描述边坡实例选取国内某矿区[13]，该边坡考虑弹性和弹塑性两种材料，边坡高约378m，坡度约为42.32º，实测经验表明，边坡的影响范围在2倍坡高范围，因此几何模型划定区域为边坡横向延伸2倍坡高，纵向延伸3倍坡高。

边坡尺寸如图2所示，边坡围岩材料属性见表2。

图2 边坡几何尺寸Fig.2 Slope geometry表2 边坡模型围岩参数Table 2 The model parameters of the rock slop类别弹性模量(Gpa)泊松比容重(KN/m3)内聚力(Mpa)内摩擦角(º)围岩2300.25250.942围岩1310.2427--3.2 计算方案及计算结果分析本次计算考虑采用双层模型，模型上部为理想弹塑性材料，模型下部为弹性材料。

这样基本消除了由于边界效应在边坡下部出现的塑性区，可以更好的模拟边坡的变形和图3 边坡有限元计算模型Fig.3 Finite element calculation mode of slope表3 强度折减参数Table 3 Reduction parameters of strength强度折减系数Fs折减后的内聚力C(Mpa)折减后的内摩擦角φ()1.00.90042.001.20.75036.881.40.64332.741.60.56329.371.80.50026.572.00.45024.232.20.40922.252.40.37520.562.60.34619.102.80.32117.823.00.30016.70塑性区的发展。

由于实际边坡纵向很长，计算模型可以简化为平面应变问题，即假定边坡所承受的外力不随Z轴变化，位移和应变都发生在自身平面内。

计算过程中边坡两侧边界水平位移受约束，下侧边界水平位移和竖向位移都受约束(见图3)，采用非关联流动法则进行计算，两层模型均选用Plane82单元，计算模型（见图3)共划分557个单元，1782个节点。

本算例中取了折减系数为1.0、1.2、1.4、1.6、1.8、2.0、2.2、2.4、2.6、2.8、3.0共11种情况对围岩2进行折减。

按式(2)进行折减，围岩2折减后的参数列入表3。

选用静力法求解，设定最大子步，打开时间步长预测器，设定牛顿-拉普森选项内容，打开大位移，线性收索设定，指定力与位移收敛精度，即可开始问题的求解[13]。

表4 求解结果Fig.4 The calculation results强度折减系数F 边坡最大变形位移(mm)边坡水平方向最大位移(mm)最大塑性应变1829.09459.31801.2830.48659.75201.4830.48659.7490.937E-061.6878.46671.3660.898E-041.8898.27772.5090.114E-032.0951.69773.9780.194E-032.21002.5371.8280.384E-032.41051.267.0940.675E-032.61096.0958.6150.862E-032.81139.3446.030.0014113.01018.3231.6020.003349求解结果表明，当折减系数为1.0时，计算机自动计算自收敛，无塑性应变(见表4）。