统一弹塑性本构模型在ABAQUS中的开发与应用

2010年4月 Rock and Soil Mechanics Apr. 2010
收稿日期：2009-05-27
第一作者简介：潘晓明，男，1979年生，博士研究生，主要从事隧道及地下结构方面。

E-mail: pxm155138@
文章编号：1000－7598 (2010) 04－1092－07
统一弹塑性本构模型在ABAQUS 中的开发与应用
潘晓明1, 2 ，孔 娟1, 2， 杨 钊1, 2， 刘 成1, 2
（1.同济大学 地下建筑与工程系，上海 200092；2.同济大学 岩土及地下工程教育部重点实验室，上海 200092）
摘 要：基于统一弹塑性本构模型的有限元理论格式，根据ABAQUS 的UMAT 格式要求，编制相应的接口程序，将统一弹塑性本构模型引入ABAQUS 中。

采用退化的统一强度模型（0b =时，为Mohr-Coulomb 模型）与ABAQUS 自带的Mohr-Coulomb 模型，对单轴试验和圆形硐室进行弹塑性分析，验证所开发材料子程序的正确性及高效性。

考虑到统一强度模型的一般情况（0b ≠）和屈服面硬化条件，对圆形硐室进行弹塑性计算，得到应力场的变化规律。

所提出的研究思路具有普遍性，为采用ABAQUS 平台进行本构模型的二次开发提供了借鉴和参考。

关 键 词：统一强度理论；屈服面；流动矢量；奇异点；应力拉回算法；ABAQUS 中图分类号：TD 353.6 文献标识码：A
Secondary development and application of unified elastoplastic
constitutive model to ABAQUS
PAN Xiao-ming 1, 2, KONG Juan 1, 2, YANG Zhao 1, 2, LIU Cheng 1,2
(1.Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China;
2.Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai 200092, China)
Abstract: Based on finite element theoretical scheme of unified elastoplastic constitutive model, and according to the UMAT interface requirement of ABAQUS, the corresponding UMAT codes are programmed, which will be called by the main analytical module of ABAQUS. Adopting degenerative model of the unified strength (0b =，Mohr-Coulomb model) and the built-in Mohr-Coulomb model of ABAQUS, the uniaxial tests and circular chamber are analyzed to verify the correctness and efficiency of the developed material subroutine. Finally, considering the general situation form of unified elastoplastic constitutive model (0b ≠) and hard condition of yield surface, which are not available in ABAQUS software, circular chamber is simulated and variational discipline of stress field is obtained. The provided basic procedures and programming essentials of the UMAT redefining in ABAQUS are universal and can offer a reference for other developers.
Key words: unified strength theory; yield surface; flow vector; singular points; return stress algorithm; ABAQUS
1 引 言
我国力学专家俞茂宏教授从多滑移单元体力学模型出发，考虑了作用在双剪应力单元体上的所有应力分量对材料屈服或破坏的不同影响，提出了一个能够适用于各种岩土类材料的统一强度理论和统一形式的数学表达式。

Mohr-Coulomb 强度理论和双剪强度理论均为其特例，并且还包含了可以比D-P 准则更合理的新的计算准则，以及可以描述非凸极限面试验结果的新的非凸强度理论[1－5]。

同时
统一强度理论可以很好地考虑中间主应力对岩体强
度的影响
[6－7]
，并能与岩土材料的真三轴试验结果
相吻合[8]，从而在岩土工程领域得到广泛的应用。

通用有限元软件ABAQUS 可以分析各种固体力学、结构力学，特别是能够处理材料非线性、几何非线性、接触非线性以及这三种非线性各种组合的高度非线性问题。

ABAQUS 本身带有较多材料非线性弹塑性模型，包括D-P 模型、Mohr-Coulomb 模型、剑桥模型等，但岩土工程领域中采用的统一弹塑性本构模型并没有在ABAQUS 中得以实现，这使得ABAQUS 在岩土工程数值分析中难以充分发挥，加入了统一弹塑性本构模型，则可以充分发挥ABAQUS 软件的计算能力。

为了弥补这一不足，本文利用用户材料子程序（UMAT ）接口[9]，通过
Fortran 或者VC 编程来开发了统一弹塑性本构模型，对单轴压缩试验的数值模拟及圆形硐室问题进行了弹塑性分析，并与ABAQUS 软件自带模型及解析解进行比较，以验证编制的接口程序的正确性。

2 统一弹塑性本构模型
将统一强度理论表示成应力不变量的形式为[10]：
(0)
b b F F θθθθθθπ=+′=
+≤≤≤≤（1）
1cos 3
θ−= （2）
tan b θ−= （3）
02cos 1sin ,1sin 1sin t c ϕϕ
ασϕϕ
−=
=++ （4）
式中：c 0为内黏聚力；ϕ为内摩擦角；b 为中间主应力系数。

如图1所示，统一强度理论包含了多个强度理论：当01b ≤≤时为统一强度理论（外凸理论）；当0b <或1b >时为双剪非凸强度理论（非凸理论）；当0b =时统一强度理论退化为Mohr-Coulomb 强度理论。

本文研究仅考虑01b ≤≤的外凸理论。

图1 统一强度理论的极限面
Fig.1 Limit surface of unified strength theory
由屈服函数()F F ′和塑性势函数()Q Q ′决定的弹塑性刚度矩阵可以写成：
[][][][][]Q F F Q A ∂∂⎡⎤⎡⎤
⎢⎥⎢⎥∂∂⎣⎦⎣⎦=−∂∂⎡⎤⎡⎤+⎢⎥⎢⎥
∂∂⎣⎦⎣⎦
T
e e ep e T
e D D σσD D D σσ （5） p
d 'd A H σ
ε== （6） 式中：A 与H ′为硬化函数；σ为等效应力；p ε为等效塑性应变。

()F F ′是统一强度理论的表达式。

如果采用()Q F Q F ′′==，由式（5）可以得到相关联流动的统一强度弹塑性刚度矩阵。

文献[10]给出了相关联流动的统一弹塑性流动矢量以及本构模型中奇异点的数学处理。

定义统一强度理论的流动矢量为
{}{}
F
∂=
∂a σ （7） {}{}
{}{}{}11123123
I F I C C C ∂∂=
+=
∂∂++a σ a a a （8）
其中：
{}{}
{}
{
}}{}{}112232
322221,1,1,0,0,0,,,2,2,2(3(),3
(),2(),3
2(x y z yz zx
xy x z yz x
z xz x y xz xz xy x yz xy yz y x I J J J J σσστττσστσστσστττστττστ∂=
=∂′′′′′′==′∂⎧′′′==−+
⎨∂⎩′′′′−+′′′′′′′′−+−′′′′−T
T a σa a σ }),2()z yz xz z xy ττστ⎫
⎪⎪⎪⎪⎪
⎪
⎪⎪⎬⎪
⎪⎪⎪⎪⎪⎪′′′′−⎪⎭
T
（9）
112321(1)3
(12sin 32F C I C C J αθθαθθ∂⎫==−⎪∂⎪⎪=⎪⎪
⎪+⎪
⎪⎪
⎬
⎤⎪
⎥⎪
⎣⎦⎪
⎪=
=⎪⎪
⎪
⎡⎤⎪
+⎢⎥⎪⎣⎦⎭
（10）
σ1
σ2
σ3
θ θb
同理，定义屈服面F ′流动矢量为 {}{}{}{}{}123123F C C C σ′
∂′′′=
=++∂a a a a （11）
123C C C θ⎫
′=⎪⎪⎪′+⎪⎪⎪⎪⎪
⎬⎪⎪⎪⎪′=
−⎪⎪⎪⎪⎦
⎭
（12） 如图1所示，统一弹塑性本构模型的屈服面应具有单一的硬化参数和流动方向，但在A 、B 、C 点存在奇异性，定义A 、B 、C 点为奇异点。

奇异点是指流动矢量{}a 在该点不能唯一确定的点。

这里定义两类奇异点，针对不同的奇异点，采用不同的处理方法。

（1）对于在b θθ=，如图1中B 点产生的奇异性，采用矢量平均的办法，即
当1tan 12b θθα−==+
1
()2
b
i i i C C C θ′=+，（1,2,3i =） （13）
（2）当1b =时，在点0θ=°和60θ=°处，如图1中A 和C 点产生的奇异性，采用数学极限的方法，确定流动矢量。

当0θ=°时，
010
20
32C C C ⎫
===⎭D D D
（14）
当60θ=°时，
6016026032C C C ′=′=′=⎭
D D D
（15） 当1b ≠时，在点0θ=°和60θ=°处产生的奇异，采用物理的方法，确定流动矢量。

当0θ=°时，
010
203
1(1)3
(120C C C αα⎫=−⎪
⎪=+⎪=⎪⎭
D D D
（16）
当60θ=°时，
6016026031(1)3
1()20C C C αα⎫
′=−⎪
⎪′=+′=⎪⎭
D D D
（17） 3 应力返回算法
施加荷载增量后，首先计算试探应力：
{}{}
[]{}1
e r r e r
−=+∆σσD ε （18）
将上述试探应力带入屈服条件（1），如果它不满足，表明此时材料的行为是弹性的。

硬化参数A 保持不变，第r 个增量步的应力就等于试探应力。

如果满足屈服条件，则表示在当前荷载步内，积分点处达到塑性条件，因此，应按照塑性规律进行计算，即该点的应力应在屈服面上移动。

分以下两种情况[11
－12]
：
第1种，如图2所示，应力状态由A 点穿越屈服面达到B 点。

此时
{}{}(,)0(,)0A B F F κκ<⎫⎪
⎬>⎪⎭
σσ （19）
图2 弹-塑性连续体增量应力变化
Fig.2 Changes of incremental stress of
elastoplastic continuum
为保持应力状态在屈服面上，则C 点的位置，可由下式决定：
A (σij ) ,F 0(F ′0)<0
A’(σij ) ,F 0(F ′0)<0
B’(σij ),F 1(F ′1)>0
B (σij ),F 1(F ′1)>0
F (F ′)=0
{}{}{}(1)e R =−−∆C A σσσ （20）
10
F R F F =
− （21） 式中：R 为比例修正因子
第2种，如图2所示，应力状态由屈服面A ′点达到B ′点。

此时
{}{}(,)0(,)0F F κκ′′=⎫⎪
⎬>⎪⎭
A B σσ （22）
1.0R = （23）
如果荷载增量比较大，且应力点位于屈服面的大曲率附近，则上述过程仍存在较大的误差。

为了提高计算的精度，可以根据超出屈服面的应力大小动态设定n 个等分数，依次将迭代应力拉回到屈服面，最后再采用上述的比例修正法进行修正，就可得到较为准确的结果。

如图3所示，将超过屈服面应力点分成n 等分，则经过n 个循环后，应力点将返回偏离屈服面的E 点，然后再通过比例修正使其
返回屈服面的E ′点。

图3 应力点返回屈服面的过程
Fig.3 Process of stress returning to yield surface
4 统一模型在ABAQUS 中的计算流程
ABAQUS 的用户材料子程序（UMAT ）通过与其求解器 Standard 的接口来实现数据交流。

UMAT 有自己的书写格式与一些规范，与主程序共享的变量必须在子程序开头予以定义。

而主程序通过ABAQUS 输入文件（.inp ）中的关键字“USER MATERIAL ”来判断用户是否使用了自定义材料本构模型，从而扩展了它的适用性和应用空间。

统一弹塑性本构模型主要求解过程[13]：每一个增量加载步开始时，ABAQUS 主程序在单元的积分点上调用UMAT 子程序，并传入应变增量、时间步长及荷载增量，同时也传入当前已知状态的应力、应变及其他与求解过程相关的变量；UMAT 子程序根据本构方程求解应力增量并更新应力及其他相关
的量，提供Jacobian 矩阵给ABAQUS 主程序以形成整体刚度矩阵；主程序结合当前荷载增量求解位移增量，继而进行平衡校核；如果不满足指定的误差，ABAQUS 将进行迭代直到认为收敛，然后进行下一增量步的求解。

根据上述思路，本文的UMAT 子程序主要流程如图4所示。

5 计算格式与程序正确性的验证
5.1 统一强度理论退化为Mohr-Coulomb 准则的
算例比较
验证算例为一个四边形平面应变单元，试件长、宽均为100 mm ，如图5所示。

首先，采用ABAQUS 自带的Mohr-Coulomb 材料，进行轴向位移压缩试验[14]；然后调用用户自定义的统一强度理论子程序，令参数b =0，使其退化为Mohr-Coulomb 强度屈服准则，同样进行单轴压缩试验；最后将两者的计算结果进行对比。

ABAQUS 软件中，应力受压为负，受拉为正。

材料参数见表1。

表1材料计算参数
Table 1 Mechanical parameters
弹性模量E /MPa
泊松比 υ
内黏聚力c /MPa
内摩擦角φ/(°)
中间主应力 系数b 200 0.3 1.0 30
图4 弹塑性分析UMAT 流程图
Fig.4 UMAT flow chart of elastoplastic analysis
图5 试件单元网格 Fig.5 Mesh of element specimen
图6为整个计算过程中单元应力σ2与轴向位移的关系曲线。

由图可见，两者计算结果十分吻合。

表明在ABAQUS 中开发实现统一强度理论（b=0）退化为Mohr-Coulomb 的计算结果是正确的。

图6 σ2与轴向位移的关系曲线
Fig.6 Relationship of σ2and axial displacement
对于统一强度理论，当b 分别取0、0.5、1.0时，单轴试验的单元应力σ2与轴向位移的关系曲线如图7所示。

从图中看出，屈服应力随b 的增加而增加，这恰与图1相吻合。

当考虑硬化参数A=1.0E 7
为常数，b 分别取0、
0.5、1.0时，单轴试验的单元应力σ2与轴向位移的关系曲线如图8所示。

单元的σ2超过初始屈服应力后，随变形的增加线性增加。

图7 不同b 值下σ2与试件轴向位移的关系曲线 Fig.7 Relationship between σ2and axial displacement
under different values of b
图8 硬化条件下，不同b 值σ2与试件
轴向位移的关系曲线
Fig.8 Relationship between σ2 and axial displacement under different values of value after hardening conditions
5.2 地下圆形硐室统一本构模型的数值解与解析
解对比
图9为一圆形硐室的计算网格，隧道半径为
10 m ，地应力p 0在x 和y 方向均为30 MPa 。

图9 硐室计算网格
Fig.9 Numerical mesh of the cave
本文取1/4模型进行计算。

统一弹塑性本构模型的计算参数如表2所示。

表2 材料计算参数
Table 2 Mechanical parameters
弹性模量E /GPa
泊松比 υ
内黏聚力c /MPa
内摩擦角φ/(°)
中间主应力 系数b
2.0
0.3 6.0 30 0,0.5,1
理想弹塑性Mohr-Coulomb 材料围岩内应力的解析表达式[15]、塑性区的径向应力p r σ和环向应力
p θσ分别为
2sin 1sin 2sin 1sin cot ()cot 1sin cot (
)(cot 1sin r r
c c a
r
c c a
ϕ
ϕϕϕ
θσϕϕ
ϕσϕϕ
ϕ−−=−−=−+ （24）
式中：a 为圆形硐室半径；r 为围岩内一点至圆心
0.00 0.01 0.020.03 0.04
-4.0
-3.5-3.0-2.5-2.0-1.5-1.0-0.50.0
解轴向位移/m
σ2/M P a
0.00
0.01 0.02
0.03 0.04
σ2/M P a
轴向位移/m
0.00
0.01
0.020.03
0.04
轴向位移/m
σ2/M P a
的半径。

设1r σ为弹塑性交界面1r r =处的径向应力，则弹性区域的径向应力e r σ和环向应力e θσ分别为
1
1
221102
2
2211022
(1)(1)e r r
e
r
r
r
p r r r r p r
r θσσσσ=−+=+− （25） 塑性区半径1r 及弹塑性交界面上的径向应力
1
r σ为
1sin 2sin 102sin 1sin 11(cot )(1sin )cot cot ()cot r p c r r c r
c c a
ϕϕ
ϕ
ϕϕϕσϕϕ
−−⎡⎤+−=⎢⎣⎦=− （26）
当0b =时，统一强度理论退化为Mohr
-Coulomb 强度理论。

如图10所示，计算所得径向应力及环向应力的数值解与解析解吻合很好。

径向应力的最大相对误差为1.55%，环向应力的最大相对误差为 1.46%，塑性区半径最大相对误差为
1.38%。

图10 径向应力及环向应力的数值解与解析解
Fig.10 Numerical result and analytical result of radial
stress and perimeter stress
图11、12为统一强度理论当0b =、0.5、1.0时的围岩内部径向应力与环向应力分布曲线。

图11 不同b 值径向应力的数值解
Fig.11 Numerical result of radial stress under
different b value
图12 不同b 值环向应力的数值解
Fig.12 Numerical result of perimeter stress under
different values
从图中可以看出，塑性区半径外，径向应力和
r σ环向应力θσ随b 值的增加而增加，但相差不大，
且随着半径的增加，最终趋于稳定值。

塑性区半径内，环向应力明显增加。

从图12可以看出，塑性区半径1r 随b 值的增加而减小，当0b =时，1 4.59r =；当1b =时，1 2.29r =。

当b=0和1时，考虑硬化参数A=5e 8为常数与不考虑硬化的径向应力及环向应力的数值解比较分别见图13、14所示。

考虑硬化后，塑性区半径有所减小，塑性区内径向应力相差不大，环向应力提高
较大。

图13 A=5e 8
, b=0径向应力与环向应力的数值解
Fig.13 Numerical result of radial stress and perimeter
stress under A=5e 8,b=0
图14 A=5e 8, b=1径向应力与环向应力的数值解 Fig.14 Numerical result of radial stress and perimeter
stress under A=5e 8,b=1
020 40 60 80 100-50
0环向应力和径向应力/M P a
沿半径方向至硐室中心的距离/m
20 40 60 80 100
0径向应力/M P a
沿半径方向至硐室中心的距离/m
20
40 60 80 100
沿半径方向至硐室中心的距离/m
环向应力/M P a
2040 60 80 100
0沿半径方向至硐室中心的距离/m
环向应力和径向应力/M P a
20
40 60 80 100
-50-40-30-20-100
沿半径方向至硐室中心的距离/m
环向应力和径向应力/M P a
6 结语
岩土材料是一种复杂的工程介质，采用现有通用有限元软件ABAQUS所提供的本构模型往往不能满足工程实际数值分析的需要。

为此，本文依据ABAQUS所提供的二次开发程序接口，结合统一弹塑性本构模型，推导了其在ABAQ-US中的增量迭代格式，并编制了相应的接口程序，可在ABAQUS 软件中进行应用。

首先通过使计算参数0
b=，将统一弹塑性本构模型退化为Mohr-Coulomb模型，进行了单轴试验模拟，并与软件自带的模型进行了对比。

结果表明，所开发的统一模型得到的结果与软件自带模型的结果十分吻合。

然后，对圆形硐室进行弹塑性分析，其计算结果与解析解吻合很好，计算精度令人满意。

最后，取b为其他值（ABAQUS所没有的本构模型），考虑强化和不考虑强化的参数，对单轴试验以及圆形硐室进行了弹塑性分析，得到了应力与位移、应力场的变化规律曲线。

以上结果表明，本文所开发计算程序是正确的，所提出的二次开发研究思路具有普遍性，对于编制自定义本构模型的程序具有一定的参考价值。

参考文献
[1]俞茂宏. 岩土类材料的统一强度理论及其应用[J]. 岩
土工程学报, 1994, 16(2): 1－9.
YU Mao-hong. Unified strength theory for geomaterials
and its applications[J]. Chinese Journal of Geotechnical
Engineering, 1994, 16(2): 1－9.
[2]俞茂宏, 昝月稳, 范文, 等. 20世纪岩石强度理论的发
展——纪念Mohr-Coulomb强度理论100周年[J]. 岩石
力学与工程学报, 2000, 19(5): 545－550.
YU Mao-hong, ZAN Yue-wen, FAN Wen, et al. Advances
in strength theory of rock in 20 century—100 years in
memory of the Mohr-Coulomb strength theory[J].
Chinese Journal of Rock Mechanics and Engineering,
2000, 19(5): 545－550.
[3]俞茂宏, 张永强, 杨松岩. 统一强度理论的又一重要推
广[J]. 西安交通大学学报, 1998, 32(12): 108－110.
YU Mao-hong, ZHANG Yong-qiang, YANG Song-yan.
Another important generalization of the unified strength
theory[J]. Journal of Xi’an Jiaotong University, 1998,
32(12): 108－110.
[4]俞茂宏, 孟晓明. 双剪弹塑性模型及其在土工问题中
的应用[J]. 岩土工程学报, 1992, 14(3): 71－75.
YU Mao-hong, MEN Xiao-ming. Application of Double-shear elastic-plastic model in the geotechnical
problems[J]. Chinese Journal of Geotechnical Engineering, 1992, 14(3): 71－75.
[5]俞茂宏. 双剪理论及其应用[M]. 北京: 科学出版社,
1998, 6: 837－843.
[6]徐德欣. 岩石中间主应力效应的理论分析[J]. 山地学
报, 2003, 21(2): 246－251.
XU De-xin. The theoretical analysis of the intermediate
principal stress of rock[J]. Journal of Mountain Science,
2003, 21(2): 246－251.
[7]李小春, 许东俊. 中间主应力对岩石强度的影响程度
和规律[J]. 岩土力学, 1991, 12(1): 9－16.
LI Xiao-chun, XU Dong-jun. Law and degree of effect
the intermediate principle stress on strength of rock[J].
Rock and Soil Mechanics, 1991, 12(1): 9－16.
[8]张金涛, 林天健. 三轴实验中岩石的应力状态和破坏
性质[J]. 力学学报, 1979, 2: 99－105.
ZHANG Jin-tao, LIN Tian-jian. Stress conditions and
variations of rupture characteristics of a rock as shown by
triaxial tests[J]. Acta Mechanica Sinica, 1979, 2: 99－
105.
[9]HIBBIT, KARLSON, SORRENSON. ABAQUS user’ s
manual[M]. Rhode Island: Pawtucket R I, 2002 .
[10]俞茂宏, 杨松岩, 范寿昌, 等. 双剪统一弹塑性本构模
型及其工程应用[J]. 岩土工程学报, 1997, 19(6): 2－10.
YU Mao-hong, YANG Song-yan, FAN Sau-cheong, et al.
Twin shear unified elastoplastic Constitutive model and
its applications[J]. Chinese Journal of Geotechnical
Engineering, 1997, 19(6): 2－10.
[11]OWEN DRJ, HINTON E. Finite elements in plasticity,
theory and practice[M]. Swansea: Pineridge Press, 1980. [12]朱伯芳. 有限单元法原理与应用[M]. 北京: 中国水利
水电出版社, 2004.
[13]杨曼娟. ABAQUS用户材料子程序开发及应用[D]. 武
汉: 华中科技大学, 2005.
[14]张传庆, 周辉, 冯夏庭. 统一弹塑性本构模型在
FLAC3D中的计算格式[J]. 岩土力学, 2008, 29(3): 596
－602.
ZHANG Chuan-qing, ZHOU Hui, FENG Xia-ting.
Numerical format of elastoplastic constitutive model
basedon the unified strength theory in FLAC3D[J]. Rock
and Soil Mechanics, 2008, 29(3): 596－602.
[15]于学馥, 郑颖人, 刘怀恒, 等. 地下工程围岩稳定分
[M]. 北京: 煤炭工业出版社, 1983, 12: 122－155.。