基于双剪统一屈服准则在岩质边坡稳定分析中的应用

Application of Double Shear Unified Yield Criterion to Stability Analysis of Rock Slope

摘要: 岩质边坡分析中常用的Drucker-Prager屈服准则，并未考虑中间主应力效应，无法反映岩石材料的真实受力状态。而双剪统一强度理论包含了一系线性破坏准则，通过改变中间主应力系数b来反映中间主应力的影响程度。利用ANSYS有限元软件基于双剪统一屈服准则与统一弹塑性有限元程序UEPP，分别计算b为0、0.5、1共三种情况下岩质边坡稳定性情况。结果表明：随着中间主应力系数b的减小，等效塑性应变分布区域增大，安全系数减小，剪应力增大，岩体的承载能力没有得到更大的发挥，不利于对实际工程优化设计。另外为了更加清晰地了解中间主应力对边坡稳定分析的影响，与Drucker-Prager屈服准则模拟结果相比较，结果表明：中间主应力在边坡稳定分析中不可忽视，并且基于双剪统一屈服准则下的岩质边坡更难贯通，能够充分发挥材料的强度潜能。

Abstract: The Drucker-Prager yield criterion commonly used in rock slope analysis does not consider the intermediate principal stress effect and cannot reflect the real stress state of rock material. The double-shear unified strength theory contains a linear failure criterion which reflects the influence degree of intermediate principal stress by changing the intermediate principal stress coefficient. Using ANSYS finite element software based on double shear unified yield criterion and unified elastic-plastic finite element program UEPP, the rock slope stability was calculated as 0, 0.5 and 1 respectively. The results show that with the decrease of intermediate principal stress coefficient, the distribution area of equivalent plastic strain increases, the shear stress increases, and the bearing capacity of rock mass does not get greater play, which is not conducive to the practical engineering optimization design. In addition, in order to understand the influence of intermediate principal stress on slope stability analysis more clearly, the simulation results are compared with Drucker-Prager yield criterion. The results show that the intermediate principal stress cannot be ignored in slope stability analysis, and the rock slope based on the unified yield criterion of double shear is more difficult to penetrate, which can make full use of the strength potential of materials.

1. 引言

边坡稳定性问题一直是矿山工程、水利水电工程、建筑工程等工程领域的难题，近几十年来，我国大力发展基础建设，作为三大自然灾害之一的滑坡灾害严重的危害我国人民的生命安全与国家财产。评价边坡稳定的方法是对解决边坡稳定问题重要的方面之一。

目前，岩质边坡稳定性分析主要是方法有极限平衡方法和有限元法，其中有限元法考虑了边坡岩体的不连续性与非均质性，避免了极限平衡方法中将滑体视为刚体过于简化的问题。大型通用有限元软件ANSYS在工业应用领域和科学研究方面都有深入的应用。其强大的分析功能能够解决各类问题，如结构、电磁、声场、热力学等 [1]。目前已被广泛应用于土木工程、水利水电、航天工程、汽车工程等近二十多个领域。但是随着在各领域的广泛发展，ANSYS也遇到了一些问题，如材料本构模型不足的问题。

岩土类材料本构关系十分复杂，其变形主要是由颗粒间的位移变化引起，不同应力水平下相同的应力增量引起的应变增量并不相同，从而表现出复杂的非线性特性 [2]。ANSYS中自带的岩土类弹塑性本构模型只有Drucker-Prager准则，它不仅忽略了岩土类材料的非线性特性，同时并没有考虑到中间主应力效应。目前，国内对ANSYS有限元软件进行材料本构模型二次主要采用的是ANSYS参数化设计语言(APDL)、用户可编程特性(UPFs)。宿辉 [3] 等采用二次开发工具APDL实现了邓肯张非线性材料本构模型的二次开发。黄强 [4] 等采用ANSYS二次开发工具APDL和UIDL，对等效DP-MC准则的有限元强度折减法进行二次开发。袁野 [5] 等采用二次开发工具UPFs对邓肯张模型进行了二次开发。关云飞 [6] 等利用二次开发工具UPFs对修正的剑桥模型进行了二次开发。众多的弹塑性有限元程序允许添加新的材料本构，如上面提到的ANSYS有限元软件 [7]。同时一些有限元程序的源程序内也可以添加新的弹塑性屈服准则，例如有限元程序UEPP [8]。

笔者将采用俞茂宏教授提出的双剪统一强度理论 [9] 与统一弹塑性有限元程序UEPP。将双剪统一屈服准则应用于实际工程的岩质边坡中，改变中间主应力系数，并将结果与Drucker-Prager屈服准则模拟结果相比较，通过分析可知，中间主应力在边坡稳定分析中不可忽视，采用双剪统一屈服准则，更能反映岩土强度特性，有效地发挥材料潜能，节约实际工程费用。

2. 基于统一强度理论的二次开发

2.1. 相关流动的双剪统一弹塑性刚度矩阵[D_ep]

1985年，俞茂宏教授提出了新的双剪单元模型，由此给出了双剪理论 [10]。在1990年俞茂宏教授增加了新的材料参数反映中间主应力效应，提出了统一强度理论，建立了一个新的双剪单元，认为作用于双剪单元上的两个较大剪应力及其面上的正应力所组成的函数到达一定的极限值时，材料发生破坏，其数学表达式为：

$\begin{array}{l}F={\tau }_{13}+b{\tau }_{12}+\beta \left({\sigma }_{13}+b{\sigma }_{12}\right)=c\\ {\tau }_{12}+\beta {\sigma }_{12}\ge {\tau }_{23}+\beta {\sigma }_{23}\\ F^{\prime }={\tau }_{13}+b{\tau }_{23}+\beta \left({\sigma }_{13}+b{\sigma }_{23}\right)=c\\ {\tau }_{12}+\beta {\sigma }_{12}\le {\tau }_{23}+\beta {\sigma }_{23}\end{array}$ (1)

式中：b为中间主应力系数；c、 $\beta$ 为材料常数。统一强度理论的 $\pi$ 平面极限线如图1所示。同时为了方便数值计算，采用Nayak应力不变量形式屈服准则，双剪统一强度理论应力不变量形式的表达式为 [11]：

Figure 1. Limit lines of unified strength theory on $\pi$ plane

图1. 统一强度理论 $\pi$ 平面图

$\begin{array}{ll}F=\left(1-\alpha \right)\frac{I_1}{3}+\frac{\alpha \left(1-b\right)}{1+b}\sqrt{J_2}\sin \theta +\left(2+\alpha \right)\sqrt{\frac{J_2}{3}}\cos \theta ={\sigma }_t,\hfill & \left(0˚\le \theta \le {\theta }_b\right)\hfill \\ F^{\prime }=\left(1-a\right)\frac{I_1}{3}+\left(a+\frac{b}{1+b}\right)\sqrt{J_2}\sin \theta +\left(\frac{2-b}{1+b}+\alpha \right)\sqrt{\frac{J_2}{3}}\cos \theta ={\sigma }_t,\hfill & \left({\theta }_b\le \theta \le 60˚\right)\hfill \end{array}$ (2)

流动矢量 $\left\{\alpha \right\}$ 各分量系数如下 [8]：

$C_1=\frac{\partial F}{\partial I_1}$ (3)

$C_2=\frac{\partial F}{\partial \sqrt{J_2}}+\frac{c\tan 3\theta }{\sqrt{J_2}}\frac{\partial F}{\partial \theta }$ (4)

$C_3=-\frac{\sqrt{3}}{2\sqrt{J_2^3}\sin 3\theta }\frac{\partial F}{\partial \theta }$ (5)

使用流动矢量 $\left\{\alpha \right\}$ 各分量系数，求解流动矢量，从而得到双剪统一弹塑性刚度矩阵 $\left[D_{ep}\right]$ 。由于双剪统一强度理论为分段线性屈服函数，当 $F=F^{\prime }$ 处会出现角点奇异性问题。基于关联流动法则，采用矢量平均的办法 [12]，来确定流动矢量 $\left\{\alpha \right\}$ ：

$F\ge F^{\prime }$ ，即 $0\le \theta \le {\theta }_0$

$\begin{array}{l}\underset{\_}{C_1}=\frac{1}{3}\left(1-\alpha \right)\\ C_2=\left(1+\frac{\alpha }{2}\right)\frac{2}{\sqrt{3}}\cos \theta +\frac{\alpha \left(1-b\right)}{1+b}\sin \theta +\cot 3\theta \left[-\left(1+\frac{\alpha }{2}\right)\frac{2}{\sqrt{3}}\sin \theta +\frac{\alpha \left(1-b\right)}{1+b}\cos \theta \right]\\ C_3=-\frac{\sqrt{3}}{2J_2\sin 3\theta }\left[-\left(1+\frac{\alpha }{2}\right)\frac{2}{\sqrt{3}}\sin \theta +\frac{\alpha \left(1-b\right)}{1+b}\cos \theta \right]\end{array}$ (6)

$F<F^{\prime }$ 即 ${\theta }_0<\theta \le \frac{\pi }{3}$

$\begin{array}{l}{C^{\prime }}_1=\frac{1}{3}\left(1-\alpha \right)\\ {C^{\prime }}_2=\left(\frac{2-b}{1+b}+\alpha \right)\frac{1}{\sqrt{3}}\cos \theta +\left(\alpha +\frac{b}{1+b}\right)\sin \theta +\cot 3\theta \left[-\left(\frac{2-b}{1+b}+\alpha \right)\frac{1}{\sqrt{3}}\sin \theta +\left(\alpha +\frac{b}{1+b}\cos \theta \right)\right]\\ {C^{\prime }}_3=-\frac{\sqrt{3}}{2J_2\sin 3\theta }\left[-\left(\frac{2-b}{1+b}+\alpha \right)\frac{1}{\sqrt{3}}\sin \theta +\left(\alpha +\frac{b}{1+b}\right)\cos \theta \right]\end{array}$ (7)

将双剪统一屈服准则，作为屈服条件，用于弹塑性有限元分析。根据以上本构模型编制的统一弹塑性材料模型子程序，在ANSYS有限元软件和UEPP-2D程序下运行，应用于岩质边坡稳定性分析中。UEPP可以适合于不同种类特性的材料，使用非常方便，该程序可以装入各种大型商业分析软件，例如ANSYS、ABAQUS、ADINA等。

2.2. 统一弹塑性有限元程序UEPP [8]

统一弹塑性程序UEPP是由俞茂宏教授等人开发的一个专门用于统一强度理论的弹塑性计算机程序，其中包含二维和三维两种版本。UEPP程序的主要功能为：1) 可以有效的解决平面应力、平面应变和空间轴对称问题。2) 采用了四边形单元，其单元节点数可以为4、8、9。3) 统一强度理论中的屈服函数，可以逼近和包含现有的屈服准则，并且适用于各种各向同性材料。4) 对角点采用光滑化处理，解决角点奇异性问题。UEPP采用先进的算法，可以对结构进行多种分析，例如弹性分析和弹性极限分析、弹塑性分析、地震响应分析等。

UEPP中采用的材料模型不只有岩土类材料，还可以适用于金属材料、塑料材料、混凝土材料等各种类型材料，因此，UEPP可以广泛的应用于化学工程、机械工程、铁道工程、航空工程、土木工程。UEPP程序包括UEPP-2D、UEPP-3D和UEPP-3DH。其中UEPP-2D程序可以对平面应力应变问题进行相应的弹塑性分析，它的荷载类型可分为集中荷载、体力和分布力三种荷载类型，采用波前法求解，同时可以求解10000个自由度以内的问题。本文将UEPP-2D程序添加到有限元软件ANSYS中进行后续的数值模拟，限于篇幅不再给出程序代码。图2为UEPP程序流程图。

Figure 2. Flowchart of UEPP program

图2. UEPP程序流程图

3. 岩质边坡稳定性数值分析

本文边坡实例采用国内新疆某矿岩质边坡，对于边坡这种纵向很长的实体，计算模型可以简化为平面应变问题。边坡有限元模型的边界条件是左右为水平约束，底面为固定约束。图3为边坡剖面图，共划分619单元，产生679个节点，采用PLANE183最新技术单元，边坡模型参数取值如表1。

Figure 3. Slope profile

图3. 边坡剖面图

Table 1. Model parameter value

表1. 模型参数取值

为了分析中间主应力大小对岩质边坡稳定性的影响，分别计算b为0、0.5、1共三种情况下边坡等效塑性应变云图、剪切应力云图，同时采用折减系数法求解不同中间主应力大小情况下安全系数大小。计算结果如图4~6。表2为安全系数比较表。图7为Drucker-Prager材料本构计算结果。通过分析可知，中间主应力对边坡的破坏效应影响越大，等效塑性应变减小，剪应力减小，边坡体的屈服强度增大，边坡体的安全系数越大，越不容易发生贯通效应，并且双剪统一屈服准则计算结果与Drucker-Prager准则的计算结果相比，D-P准则下的贯通性更明显，边坡稳定性偏安全。由此可见，在实际工程中，中间主应力对边坡稳定性分析不容小觑。

(a) 等效塑性应变云图 (b) XY方向切应力

Figure 4. b = 0 the equivalent plastic strain cloud and XY direction shear stress cloud

图4. b = 0时，等效塑性应变云图与XY方向切应力云图

(a) 等效塑性应变图 (b) XY方向切应力

Figure 5. b = 0.5, the equivalent plastic strain cloud and XY direction shear stress cloud

图5. b = 0.5时，等效塑性应变云图与XY方向切应力云图

(a) 等效塑性应变图 (b) XY方向切应力

Figure 6. b = 1, the equivalent plastic strain cloud and XY direction shear stress cloud

图6. b = 1时，等效塑性应变云图与XY方向切应力云图

Figure 7. Equivalent plastic strain diagram for Drucker-Prager material

图7. Drucker-Prager材料本构的等效塑性应变图

Table 2. Comparison table of safety factor results

表2. 安全系数结果比较表

4. 结论

岩质边坡稳定性分析中常用Drucker-Prager与Mohr-Coulomb屈服准则，但是它们忽略了中间主应力的影响，无法反映岩土类材料的特性。本文采用俞茂宏教授提出的统一弹塑性有限元程序(UEPP)，结合ANSYS有限元软件，将双剪统一屈服准则加入到边坡稳定性分析当中，为了分析中间主应力对边坡稳定性分析的影响程度，计算b为0、0.5、1三种情况下等效塑性应变、XY方向剪切应力大小，并且利用折减系数法求解不同主应力大小情况下安全系数大小。同时本文以Drucker-Prager准则为例，讨论不同准则对岩质边坡稳定性影响。得到如下结论：

1) 随着中间主应力对边坡的破坏影响增大，边坡体的等效塑性应变与XY方向切应力变小，屈服强度增大，安全系数增大，边坡体不易贯通破坏。

2) 将不同主应力下的计算结果与Drucker-Prager准则的计算结果相比，安全系数偏小，D-P准则下的贯通性更明显，由此表明D-P强度准则下分析得到的边坡稳定性偏安全。可见D-P本构模型并没有充分发挥材料的潜力，在实际工程中会产生巨大的浪费。

NOTES

^*通讯作者。