非线性有限元法中的屈服准则

五种屈服准则在非线性有限元法中的应用比较及问题
李珍 1 章青 1 吴旭东 1
（1．河海大学工程力学系， 南京 210098）
摘要：目前屈服准则的选用日益受到人们的重视。本文较全面地分析了适用于岩、土和混
凝土的五种屈服准则的表达式， 推导了在非线性有限元程序及商业程序 ANSYS 中的实现。 通 过分析一个简支梁应力应变规律及某重力坝沿建基面的抗滑稳定性， 对各种准则进行了较全 面的比较，得出了一系列重要规律，并发现了一些有待解决的问题,为岩土和混凝土材料非 线性问题中屈服准则的选用提供了新的依据。
不同方法得到的安全系数 D-P2 3.17 D-P3 2.91 D-P4 3.09 抗剪断强度公式 2.86 方法 安全系数 D-P1 3.74
由此看出，用 D-P3 算得的安全系数与用抗剪断强度公式得到的最为接近，误差为 1.7 ％，D-P4 相差稍大，D-P2 相差较大，而 D-P1 相差最大，误差为 30.7％，说明现在常用的 屈服准则即 D-P1 太过保守，与规范规定的刚体极限平衡法差别太大。这一规律与算例 3.1 得到的基本一致。因此在稳定性计算中可以选用 D-P3 准则，与文献[5]的结论有所不同。 另外在计算过程中，发现了一个问题：如果按照算例 3.1 的规律，由 D-P2 得到的安全 系数应该与用 D-P3 的很接近且比 D-P4 的要小，而本例中却恰好相反。这说明由 D-P 准则 的四种形式得到的安全系数并不是对所有的材料参数都满足以上规律， 而是与 c, ϕ 的范围有 关。这从各种准则的表达式中也可以看出，各种准则对应的 α , k 并不是 c, ϕ 的线性函数，因
算例 准则
M-C 0.239
D-P1 0.405
M-C D-P1 D-P2 D-P3 D-P4
D-P2 0.238
4.E-04 3.E-04 3.E-04
D-P3 0.237
D-P4 0.271
简支梁
2.2 1.7 应力值(Mpa) 1.2 0.7 0.2 0 10
位移值(m)
2.E-04 2.E-04 1.E-04 5.E-05
Dp = ∂f ∂f 其中： A = − B, D ∂σ ∂σ
T
1 ∂f ∂f 1 T D D = D aa D A ∂σ ∂σ A B=
T
（3）
∂f T ∂f ， w p 为塑性功，f 为所选的屈服函数，D 为 σ ∂σ ∂w p
弹性矩阵。此时需要一种适合于数值计算的形式来表示矢量 a = 函数[7]，因而有：
3 (3 + sin ϕ )
3 (3 − sin ϕ ) 3 (3 + sin ϕ )
内角点外接圆
内切圆
(3 + sin ϕ )
2
(3 + sin ϕ )
2
作者简介：李珍，女，1982 年生，硕士研究生。研究方向为大坝稳定及非线性数值算法。
D-P4
等面积转换圆
2 3 sin ϕ
2 3π 9 − sin 2 ϕ
关键词：抗滑稳定安全系数；非线性有限元法；弹塑性问题；屈服准则；ANSYS
对大量的工程问题，有限单元法已经成为一种强有力的数值解法。其中在土坝、岩土地 基、重力坝沿坝基的抗滑稳定性和加固、地下洞室和边坡的稳定性分析等诸多问题中，由于 【 】 岩、土和混凝土材料具有典型的材料非线性性质 1 ，所以这些分析计算都应按照材料非线 【 】 性有限元方法来解决 2 。这些问题的关键在于选用适当的屈服准则判断结构所处的状态， 进而采用不同状态下对应的本构关系进行分析计算。 1 屈服准则的种类及其表达式 工程中最为常用是摩尔－库仑屈服准则（以下简称 M-C） ：
(
)
6 3c cos ϕ
2 9 − sin 2 ϕ
(
)
2 各种屈服准则在非线性有限元程序中的实现 2.1 屈服准则在用 Fortran 语言编制的程序中的实现 由非线性有限元的基本理论和计算过程[6]知屈服准则在其中的作用有两个方面， 其一是 根据所选的屈服准则判断一点是否进入塑性变形状态。 本文通过将有效应力和等效屈服应力 相比较来判断。各种屈服准则相应的有效应力和等效屈服应力分别为式（1）和式（2）的左 边项及右边项[7]，其中 D-P1~D-P4 中的 α 、k 分别见表 1 中的公式。如果左边项的值达到右 边项的值，则必定发生塑性变形。 其二是计算与所选屈服准则对应的本构关系，并在程序中实现。本文考虑的是弹塑性 问题，因此主要任务是求塑性矩阵 D p ，其数学表达式如下[8]：
【5】
，见表（1） 。
图(a) 各屈服准则在π平面上的曲线
表 1 各种关系下用 c, ϕ 表示 α , k 的公式 编号 D-P1 D-P2 D-P3 准则特点 外角点外接圆
α
2 sin ϕ
2 sin ϕ
sin ϕ 3
k
6c cos ϕ 6c cos ϕ
3c cos ϕ 3
3 (3 − sin ϕ )
E = 2.5 × 10 4 Mpa ， 硬化参
数 H′ = 0 ， 泊 松 比
图(b) 简支梁单元划分示意图
µ = 0.3 ，内摩擦角 ϕ = 50 ° ，
粘结力 c = 2.0 Mpa 。其中长度单位为 m，应力单位为 Mpa，不考虑自重。对这个算例进行 了两个方面的工作： （1） 、计算应用不同屈服准时对应的结构所能承受的极限荷载。这里所说的“极限荷载”是 指使程序收敛的最大荷载。结果如表 4 所示： 表 4 极限荷载对比（单位：Mpa）
1 1 I 1 sin ϕ + cos θ − sin θ sin ϕ J 2 = ccos ϕ 3 3
（1）
式中： I 1 为应力张量第一不变量， J 2 为应力偏量第二不变量， θ 为应力罗代角。 式（1）在主应力空间的屈服面是不规则的六棱锥面且有一个奇异的顶点 3 ，导致数值 【 】 计算的不便甚至不收敛，因此常常采用 Drucker-Prager 提出的屈服准则 4 来进行修正，他
图 （c）重力坝网格剖分图 表 5 各部分材料参数
材料号 参数
c(凝聚力:Mp) 1.0 2.0 0.5
ϕ (内摩擦角:度)
48 55 35 表6
应变值
E(弹模:Mp) 20000 25000 25000
ν
(泊松比) 0.2 0.18 0.18
ρ (密度:kg/m3)
2500 0 0
1 2 3
其中坝体材料号为 1，地基材料号为 2，建基面处的材料号为 3，硬化参数均为零。
2 sin ϕ '
3 3 − sin ϕ ' = 2 sin ϕ
(
)
3 (3 + sin ϕ )
（7） （8）
6c ' cos ϕ '
3 3 − sin ϕ ' = 6c cos ϕ
(
)
3 (3 + sin ϕ )
将由（7） 、 （8）二式求出一组 c ' , ϕ ' 值填入 ANSYS 中的输入数据中就可以得到 D-P2 屈服准 则。依此类推得到对应于 D-P3、D-P4 的 c ' , ϕ ' ，应用十分简便。 3 算例及规律 3.1 对简支梁应力应变规律的分析对比 根据 2.1 中的讨论，编制非线性有限元程序，对如图（b）所示的简支梁进行计算分析， 梁长 4 米，高 0.6 米，顶部正中受集中荷载，强度为 p=0.23Mpa。按平面应力问题考虑，采 用平面 4 节点等参元， 单元划 分见图（b） 。弹性模量
20 30 长度(m) 图1 单元形心最大主应力
40
0.E+00 0 1 2长度(m)3 图2 节点x方向位移 4
M-C D-P1 D-P2 D-P3 D-P4
5
作者简介：李珍，女，1982 年生，硕士研究生。研究方向为大坝稳定及非线性数值算法。
9.E-05 （2） 、计算在同样的荷载作用下，应用各 M-C 8.E-05 种屈服准则所得到的应力、位移和有效塑 D-P1 7.E-05 D-P2 性应变。 取梁的底层的第 1～10 个单元 （图 6.E-05 D-P3 1～图 3）进行分析。 5.E-05 D-P4 通过表 4 和图 1～图 3，可以看出： 4.E-05 1、 D-P2 屈服准则与 M-C 准则计算得到的 3.E-05 结果最为接近，D-P3 也十分接近； 2.E-05 2、应用 D-P1 得到的位移和有效塑性应变 1.E-05 0.E+00 最小，应力值最大，D-P4 次之。说明 0 10 20 30 40 长度(m) 在同样的荷载下应用 D-P1 判断结构最 图3 单元形心有效塑性应变 后达到屈服， 屈服区最小， D-P4 次之。 3、应用 D-P3 时得到的节点位移和有效塑性应变最大，应力最小，说明应用 D-P3 时结构最 早进入塑性阶段，D-P2 结果与之十分接近。 4、对同样的结构应用 D-P1 准则得到的极限荷载最大， D-P2 较小，D-P3 最小。D-P1 比 D-P3 大了约 70％，相差很大。 5、不考虑有限元计算，仅将算例中 c, ϕ 值代入各种准则的表达式中，求解各组 α , k 值。发
作者简介：李珍，女，1982 年生，硕士研究生。研究方向为大坝稳定及非线性数值算法。
服准则，即 D-P2、D-P3、D-P4 在 ANSYS 中的实现。基本思想为：由于各种屈服准则中不 同的只是通式（2）中的参数 α , k ，而 α , k 又是材料参数 c, ϕ 的函数，恰好 ANSYS 中 D-P 准则正是通过输入 c, ϕ 值实现的， 那么可以假设 c, ϕ 为真实的材料参数， 只需通过换算求出 一组 c ' , ϕ ' 值，将其输入 ANSYS 中就能够得到想要的屈服准则类型。比如选用 D-P2 时，可 按下式进行换算：
其中：a 1 =