当前位置:首页 >> 电力/水利 >>

Hoek-Brown准则的主应力回映算法及其二次开发


第 32 卷第 7 期 2011 年 7 月
文章编号: 1000-7598 (2011) 07-2211-08

岩 土 力 学 Rock and Soil Mechanics

Vol.32 No. 7 Jul. 2011

Hoek-Brown 准则的主应力回映算法及其二次开发
陈培帅 1,陈卫忠 1,贾善坡 1, 2,杨建平 1
(1. 中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室,武汉 430071;2. 长江大学 城市建设学院,湖北 荆州 434023)

摘 要:在弹塑性模型数值实现的过程中,需要进行应力更新的回映算法。针对三维应力空间回映算法在奇异点收敛性方面 的不足,提出主应力空间的回映算法,讨论了算法实现过程的应力空间转化问题,分析了应力更新过程中确定回映区域的方 法,建立了相应的一致性刚度矩阵。基于大型有限元软件 ABAQUS 提供的用户材料子程序 UMAT 接口,利用 Fortran 编程 语言,实现了 Hoek-Brown 准则主应力空间的隐式积分算法,利用开发的模型,进行了岩石常规三轴压缩试验的数值模拟, 通过与 ABAQUS 内部的 Mohr-Coulomb 准则计算结果的对比,验证了模型和程序的可行性和准确性。 关 键 词:霍克-布朗本构模型;ABAQUS;用户材料子程序;主应力空间回映算法 中图分类号:TU 45 文献标识码:A

Stress return mapping algorithm of Hoek-Brown criterion in principal stress space and its redevelopment
CHEN Pei-shuai1,CHEN Wei-zhong1,JIA Shan-po1, 2,YANG Jian-ping1
(1. State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China; 2. School of Urban Construction, Yangtze University, Jingzhou, Hubei 434023, China)

Abstract: In the course of implementing elastoplastic model in the FEM, it is necessary to carry out the process of stress return mapping algorithm. It may be not converging at singular points in three-dimensional space. The transforming method in principal stress space is put forward to overcoming this shortcoming. Also, the process of determining the region of return mapping is analyzed. At last, the consistent stiffness matrix is deduced. Based on the implicit algorithm in principal stress space, the UMAT subroutine of the Hoek-Brown model is developed in ABAQUS with the Fortran programming language. Through the example of triaxial compression test, it is easy to verify the feasibility and accuracy of this method by comparing the result with Mohr-Coulomb criterion in ABAQUS. Key words: Hoek-Brown criterion; ABAQUS; UMAT; return mapping algorithm in principal stress space

1





综合考虑了岩块和结构面的强度以及岩体结构等多 种因素的影响,能够反映岩体的非线性破坏特征, 较好地解决了 M-C 强度准则在低应力区不太适用和 受拉破坏处理上的困难, 并已建立了一套比较完善的 参数确定方法,因而 H-B 准则更适宜于描述岩体的 破坏特征,并在岩体工程中得到了广泛的应用 。 在有限元数值模拟中,H-B 模型尚未得到广泛 的应用, 这是由于 H-B 模型存在棱线和尖点处的奇 异点,在这些奇异点处导数不连续,再加上 H-B 准 则的非线性特征,导致数值实现存在一定的难度。 [8] 为解决收敛性问题,Wan 利用圆滑函数对 H-B 准 则的奇异点进行圆滑处理,并提出了数值实现的方
[5-7]

岩土工程稳定性分析一般采用 Mohr-Coulomb 准则 (简称 M-C 准则) 或者 Drucker-Prager 准则 (简 称 D-P 准则)等屈服准则作为衡量岩土体破坏的标 准,这些准则的强度包线都是线性的,很多已有的 稳定性分析成果和经验都是建立于这些准则之上 的
[1-3]

。事实上,岩土材料的强度包线并非呈线性,

对于岩体工程而言,采用非线性的强度包线更符合 实际。 Hoek-Brown 准则(简称 H-B 准则)是 Hoek 等 人于 1980 年提出的非线性强度准则 , 由于该准则
[4]

收稿日期:2010-02-01 基金项目:国家高技术研究发展专项经费(863计划)(No. 2007AA11Z108);中国科学院“ 西部行动计划高新技术”项目( No. KGCX2-YW-506)。 第一作者简介:陈培帅,男,1985 年生,硕士研究生,主要从事隧道及地下工程施工及围岩稳定方面的工作。E-mail: chenpeishuai@foxmail.com

2212









2011 年

[9] 法。Crisfield 利用 D-P 准则对奇异点进行估量运

平衡迭代收敛的速度有直接的联系,目前主要的算 法有两种:显式算法和隐式算法。显式算法的稳定 性是有条件的, 必须在计算时恰当地控制积分步长, 否则将导致比较大的计算误差。相对而言,隐式算 法是无条件稳定的,本文采用目前常用的隐式向后 欧拉积分算法,实现计算过程中的应力更新。贾善 [11] 坡等 详细论述了算法的方法和步骤, 如图 1 所示。 其基本原理是:首先是弹性预测阶段,假设增量步 的应变增量全部为弹性应变,计算出应力更新后的 数值;其次是塑性修正阶段,检查弹性预测阶段的 应力更新数值是否满足屈服条件,若不满足屈服条 件,则无需进行修正,反之则对应力进行修正,使 应力返回到屈服面。

算,解决导数不连续问题。然而这些方法的计算结 果同经典的 H-B 模型的结果又有一定的误差, 本文 论述 H-B 模型主应力空间的隐式积分算法, 并以大 型有限元软件 ABAQUS 为平台,采用 Fortran 语言 编制了 UMAT 用户材料子程序, 通过常规三轴压缩 试验的数值模拟验证了程序的精确性和有效性。

2

H-B 破坏准则
H-B 准则是 Hoek 等通过对大量岩石三轴试验

资料和岩体现场试验成果的统计分析,并综合考虑 了岩体结构、 岩块强度、 应力状态等多方面的影响, 而得出的岩体破坏经验公式:
? ? ? f ? ? 1 ? ? 3 ? ? ci ? s-mb 3 ? ? ci ? ?
a

(1)
?n
f >0 ??
e

式中:? 1 、? 3 分别为破坏时的最大和最小有效主应 ? ci 力, 采用岩土工程中的符号惯例, 以压应力为正; 为岩石的单轴抗压强度; mb 、 a 、 s 均为经验性参 数, a 反映材料的非线性程度,当 a =1 时,H-B 准 则转化为线性强度准则;mb 为破裂的岩性参数, 与 颗粒的咬合程度有关, 随着破损程度的增加而减小;

f =0

- ?? p

?
f <0

? ??

? n+1 ?

a 、s 为与岩体特征有关的参数,s 与岩样中存在的 断裂程度有关,是岩石黏聚力的代表,对于完全无 破损的岩块,它的值是 1,对于高度破损的岩石, 它的值趋近于 0,其数值由下式确定:
? mb =mi e ? ? (GSI -100) (9-3D) s =e (2) ? ? 1 1 a = + (e ? GSI 15 ? e?20 3 ) ? 2 6 ? [10] 式中 : GSI 为 Hoek 等 提出 的地质强度 指标
(GSI -100) (28-14D)

图 1 回映算法的原理 Fig.1 The principle of return mapping algorithm

在应力更新过程中需要求屈服函数的一阶导以 及塑性势函数的一阶导和二阶导,H-B 屈服准则是 以主应力的形式描述的, 因而相对于三维应力空间, 在主应力空间实现这些求导过程更为简便
[12]

。对于

各向同性问题,应力更新和刚度矩阵的计算可以不 考虑所采用的坐标系统,这是因为在回映过程中, 剪应力始终为 0,主应力的方向不会改变,在主应 力空间求出更新值之后,根据之前的方向矩阵将应 力、 应变和一致性刚度矩阵再转化到三维应力空间。 下面详细论述主应力空间 H-B 准则隐式积分算法的 实现方法。 2.1 试探应力的求解及应力空间的转化 假定增量步开始时应力为 ? n , 当前增量步的应 变增量为 Δ? n ,设弹性应变增量等于总应变增量, 即塑性应变增量为 0,计算弹性试探应力 ? tr ,即

(geological strength index) ,其值取决于岩体结构、 岩块嵌入状态和岩体中不连续面质量等因素,也可 由岩体的 RMR 分值确定。 取塑性势函数与屈服函数的表达式一致,即
? ? ? g ? ? 1 ? ? 3 ? ? ci ? sg -mg 3 ? ? ci ? ?
ag

( 3)

当式中的参数与屈服函数中对应的参数相等 时,即 a = ag , m = mg , s = sg ,则为相关联流动; 反之则为非相关联流动。

? tr ? ? n ? DΔ? n

(4)

3

主应力空间的隐式积分算法

构造方向矩阵 H,将试探应力通过方向矩阵 H 转化到主应力空间,如下式:
?? tr 0 ? ? H ? tr ? H ?? ? ? ? ? ? ? ? x y z xy xz yz ? ?1?3 1?3? ? ?

在数值实现的过程中,应力更新是核心的环节 之一,应力更新算法的选取同计算误差的大小以及

( 5)

第7期
tr tr

陈培帅等:Hoek-Brown 准则的主应力回映算法及其二次开发

2213

将求出的 ? 1 和 ? 3 代入式(1) ,若 f ≤ 0 ,则 无需进行塑性修正,弹性试探应力即为更新的应力 值,若 f ? 0 ,则需要下一步的塑性修正,在塑性 修正的过程中,需要利用弹性刚度矩阵,由于是在 主应力空间进行回映计算,所以只需要利用刚度矩 阵中对应主应力的分量,即下式的 D :
?D 3?3 D?? ?0 ?3?3 0? 3?3 ? G? 3?3 ?

拉应力区

?3 ?1 >?2 >?3

?1 ?1 >?2 >?3 拉应力棱线 ?1 =?2

?1
压应力区

?1 >?2 >?3
压应力棱线 ?1 =?2

( 6)

2.2 确定回映区域 在研究中,我们考虑的是各向同性问题,因而 只需要研究屈服面的六分之一,研究区域的应力点 满足 ? 1 ≥ ? 2 ≥ ? 3 ,如图 2 所示的区域。

图 3 回映区域的划分 Fig.3 Division of the return mapping region

2.3 应力更新 在应力更新过程中,需要用到屈服函数的一阶 导数 f? 以及塑性势函数的一阶导数 r 和二阶导数 c ,这在主应力空间很容易求解:
?f ? ?k ?? ?g r? ? ?k g ?? ? dk g ? d? 1 2 ? g ? c? ?? 0 ?? 2 ? ? 0 ? ? f? ? 0 ? 1? 0 ? 1?
T

?3

(7) (8)

T

拉应力棱线

?1 ? =0

-30° ? 压 应 力 棱 线

?2

?1 >?2 >?3

? 0 0? ? 0 0? ? 0 0? ? ?
a ?1

(9)

图 2 垂直于静水轴的横截面在? 平面投影 Fig.2 Projection of the cross-section which is perpendicular to hydrostatic axis on plane of ?

其中:
? ? ? a g ?1 ? ? ? ? ? ?g kg ? ? 1 ? a g mg ? s g ? mg 1 ? ? (10) ?? 1 ? ci ? ? ? ag ? 2 ? 2 dk g a g mg ? ?1 ? ? ? ?1 ? a g ? ? sg ? mg ? ? d? 1 ? ci ? ? ci ? ? n ?1 如图 4 所示,取屈服面上一点 ? 作为应力回 ? ? ? ?f k? ? 1 ? am b ? s ? m b 1 ? ?? 1 ? ci ? ?

如图 3 所示,应力回映区域的有 4 种形式: ①回映到面;②回映到压应力棱线 l1 ,即 ? 1 ? ? 2 ≥

? 3 ;③回映到拉应力棱线 l2 ,即 ? 1 ≥ ? 2 ? ? 3 ; ④回应带尖点处,即 ? 1 ? ? 2 ? ? 3 。
根据弹性试探应力,从而确定回映区域分为以 下几个步骤: (1)判断是否需要回映到尖点处,如果不满足 回映到尖点的条件,则转到第 2 步; (2)直接回映到屈服面上; (3)分析第 2 步求出的应力更新值; ①若 ? 1 ? ? 2 ? ? 3 , 说明采用的回映到面的算法 是正确的,直接进行一致性刚度矩阵的求解; ②若 ? 2 ≥ ? 1 , 说明采用的回映到面上的算法是 不正确的,需要回映到压应力棱线 l1 ,即在回映迭 代的过程中始终保持 ? 1 ? ? 2 ; ③若 ? 3 ≥ ? 2 ,说明采用的回映到面上的算法 是不正确的,需要回映到拉应力棱线 l2 ,即在回映 迭代的过程中始终保持 ? 2 ? ? 3 。

映的初始值,回映的目标就是在屈服面上不断移动

? n ?1 ,使其满足应力迭代收敛的条件。按照正交化 法则, ? n ?1 处应力塑性修正方向的理论值为
??1 ? ? ? kg ? ? ? ? ? E u ? Dr ? ?? kg ? ? ? ?1 ? ? ??1 ? 2? ? ? ? ?? kg ? 1 ? ? ?

(11)

则 ? n ?1 处应力塑性修正方向在 ? 1 ? ? 3 应力空 间的斜率的理论值为
j1 ? u3 u1 ?

? kg ? 1 ? ?

?1 ? ? ? kg ? ?

(12)

在一个迭代步中, 连接 ? n ?1 和 ? n 的斜率, 即塑

2214









2011 年

性修正方向在 ? 1 ? ? 3 应力空间斜率的计算值为
j2 ?
n ?1 n ?3 ? ?3 ? 1n ?1 ? ? 1n

D con ? D c ?

D c r ? f? ? D c

T

(13)

? f? ?

T

Dc (r )

(20)

式中: D c 为修正的弹性刚度矩阵,且有:
D c ? TD
f =0

?3 ?1

(21)

其中修正矩阵为
?n
j2 1 1 j1

? in ?1
f?
g(ag=1)

u

? ?2 g ? T ? ? I ? ?? D 2 ? ?? ? ?

?1

(22)

r

上面的一致性刚度矩阵的求解公式只实用于面 回映的情况, 不适合于线回映和尖点回映。 如式 (22) 所示,求解修正矩阵 T 的过程中,需要求解塑性势 函数的二阶导,在三维应力空间中,这个求解的过 程相当复杂,而在主应力空间中,求解要相对简单 的多,如式(9)所示。 2.4.1 面回映的一致性刚度矩阵 首先计算修正矩阵 T 中对应正应力的分量:
T ? T1? 3,1?3
n ?1 1

图 4 应力更新示意图 Fig.4 Sketch of stress update

令 Q ? j2 ? j1 ,则迭代收敛的条件为: Q ? 0 , n ?1 而初始的 ? 1 一般是不满足此收敛条件的,采用 n ?1 Newton-Rapthon 方法进行迭代计算,对 Q 在 ? 1 , i 处 进行 Taylor 展开:
n ?1 Q ? Q ? 1, ? i n ?1 dQ ? 1, i n ?1 d? 1, i

?

?

?

? ??
n ?1 1, i

(23)

?0

(14)

对应的塑性乘子为
?? ? ?? p Dr

由上式可得
n ?1 n ?1 ? 1, i ?1 ? ? 1,i ?

(24)

? ? Q ? ?? ?
Q ?
n ?1 1,i

(15)

从而得到:
? ?2 g ? T ? ? I ? ?? D 2 ? ?? ? ? 3?3
?1

当满足下式的迭代收敛条件时,退出迭代步:
n ?1 n ?1 ? 1, ? TOL i ?1 ? ? 1,i

(25)

(16)

式中:TOL 为迭代收敛误差,本程序实现过程中取 TOL 为 1 ? 10?5 。 n ?1 将满足迭代条件的 ? 1 , i 代入 f ? 0 ,求解第三 主应力:
1 ? ? ? 1n,? i ? ? ci ? s ? mb ? ? ? ci ? ? ? a

然后求解修正矩阵 T 对应剪应力分量的部分 TG ? T4? 6,4? 6 ,其求解过程更为简便:
p ? ? ?? 1p ? ?? 2 ?? 0 0 ? ? n ?1 ?? n ?1 ? ? ?1 ? ? 2 ?? ? p p ? ?? ? ? ? ? ? 1 2 ?? TG ? ? I ? ? 0 0 n ?1 ? 3?3 ? ? 1n ?1 ? ? 2 ?? ? ? p p ?? ?? 1 ? ?? 2 ? ? ? ? 0 0 n ?1 ? ? ? ? ? 1n ?1 ? ? 2 ? ?? ? ?1

?

n ?1 3 ,i

??

n ?1 1 ,i

(17)

中间主应力通过下式进行求解:

? n ?1 ? ? n ? ? u

(18) 至此,可以得到修正矩阵:
?T 0 ? T=? ? ?0 T G ? ? ?

(26)

? 1n ?1 ? ? 1n 式中: ? ? ,则有: u1
n ?1 n ?2 , i ? ? u2 ? ? 2 , i

(19)

(27)

2.4 求解一致性刚度矩阵 一致性刚度矩阵不影响应力更新的最终结果, 但与整体平衡迭代的收敛速率有直接的联系,在三 维应力空间中,一致性刚度矩阵需要奇异于塑性应 变的修正方向,其求解式为

将 T 代入式 (21) 求出弹性刚度矩阵的修正值, 相应面回映的一致性刚度矩阵可以求得
Dcon ? D c ? Dc rf?T D c f?T D c r

(28)

第7期
T

陈培帅等:Hoek-Brown 准则的主应力回映算法及其二次开发
T

2215

T T 式中: f? ? ? ? f? 0 0 0 ? ? ,r ? ? ?r 0 0 0? ? 2.4.2 线回映的一致性刚度矩阵

TG 的表达方式和面回映相同,如式(26) 。从而得

到修正矩阵 T,代入式(21)求得弹性刚度矩阵的 修正值:
D c ? TD

相对于面回映而言,线回映的修正矩阵需要进 行一些修改,其对应正应力的分量为
? ? 2 gn ? ?2 g Tl ? ? I ? ?? D 2 ? ??n D 2 ? ?? ?? ? ?
?1

(38)

(29)

类似于面回映的情况,定义垂直于势函数线方 向向量的任意向量,则线回映的一致性刚度矩阵必 须奇异于这个向量,参照文献[9],一致性刚度矩阵 的表达形式为
? n ng ? ? ? ? n ?T D c ? ? 0 ?

?2 gn 式中: 2 为压应力区或者拉应力区对应的塑性势 ?? 函数的二阶导。 压应力区对应的区域满足 ? 2 ? ? 1 ? ? 3 ,则 T ?g rn ? n ? ? 0 kg ? 1? (30) ? ? ??
0? ?0 0 ? ? dk g ? gn ? ? ? 0 0 ? d? 1 ?? 2 ? ? ? ?0 0 0?
2

D con

? ? ? ?

T ?1

ng

? 0? ? ? c? G ?

(39)

(31)

c 式中:D c 和 G 分别为 D c 左上角和右下角的 3 ? 3 矩 阵;n 和 n g 分别为对应回映的棱线的屈服函数对应

拉应力区对应的区域满足 ? 1 ? ? 3 ? ? 2 ,则 T ?g rn ? n ? ? kg ? 1 0 ? (32) ? ? ??
? dk g ? 0 0? ? d? 1 ? ? 2 gn ? ? ? ? 0 0 0 2 ? ? ?? 0 0? ? 0 ? ? ? ?

的方向向量和势函数对应的方向向量。 压应力棱线满足 ? 1 ? ? 2 ,即
?? 1 ? ? ? ?? ?? 2 ? ? ? ?? 3 ? ?

(33)

?1 ? ? ? ? ?1 ? ? a? ? ?1 ? ? ?? ? ? ? ? 1 ci ? s ? mb ? ? ci ? ? ? ? ?

(40)

对应的屈服函数和势函数的方向向量分别为
n ? ?1 1 k ?
T

回映到线时,塑性应变增量的方向可以表示为 两个相邻势函数面的法向量的线性组合:
?? ? ?? r ? ??n rn
p

(41)
T

(34)

ng ? ? ?1 1 kg ? ?

(42)

当回映到压应力线时
???1p ? ? kg ? ?0? ? p? ? ? ? ? p ?? ? ??? 2 ? ? ?? ? 0 ? ? ??n ? kg ? ? p? ? ?1? ? ?1 ? ? ? ? ?? 3 ? ? ? ? ?

拉应力棱线满足 ? 2 ? ? 3 ,即

(35)

得到塑性乘子
?? ?
p ??1p ?? 2 ,??n ? kg kg

(36)

?1 ? ? ? a? ? ?1 ? ? ?? 1 ? ? ? ? ?? 1 ? ? ci ? s ? mb ? ? ? ? ? ?? 2 ? ? ? ci ? ? ? a? ? ? ? ? ? 3? ? ? ? ? ? 1 ? ? ci ? s ? mb 1 ? ? ? ? ci ? ? ? ? ?

(43)

当回映到拉应力线时,有:
???1p ? ? kg ? ? kg ? ? p? p ? ? ? ? ?ε ? ??? 2 ? ? ?? ? 0 ? ? ??n ??1 ? ? p? ??1? ?0 ? ? ?? 3 ? ? ? ? ? ? ?

对应屈服函数和势函数的方向向量分别为
n2 ? ?1 k k ?
T

(44)
T

(37)
g n2 ?? ?1 kg kg ? ?

(45)

p p 得到塑性乘子 ?? ? ??? 3 , ??n ? ??? 2 ,代入 式子 (29) 即可求得线回映修正矩阵的正应力分量。

2.4.3 点回映的一致性刚度矩阵 若应力回映到尖点处,则一致性刚度矩阵必须

2216









2011 年

垂直于任意过尖点的势函数面的法向方向,由 H-B 势函数面特征得

n2?1 = u p ? l ,取棱线上任一点 ? l ,弹性试探应力为
tr ? tr ,定义方向 (? ? ? l ) ,则有: ①若 n2?1 (? tr ? ? l ) ≤ 0 , 对应第 1 种情况, 回映 到面。 ②若 n2 ?1 (? tr ? ? l ) ? 0 ,对应第 2 种情况,回映 到棱线。

Dcon ? 0
2.5 回映结果转换到三维应力空间

(46)

利用之前的方向转化矩阵 H 将计算的更新值 转换到三维应力空间:
T

? ? H ?1 ?? tr 0 ? 3?3 ? ?
D
con xyz

(47)
H
?1

?3

?? p n2 n1

? tr

? H

?

?1 T

?

D

con

(48)

u2 ut

u1

4

尖点回映的条件
在分析尖点回映的条件之前,先引入一个边界
?1

?2

面的概念,为讨论方便,采用线性的屈服准则(棱 线是直线而非曲线) 进行分析, 如图 5 所示。 向量 u p 和 l 组成的面为 1/6 屈服面的一部分,l 为棱线,其 方向向量为 l ,取棱线上任一点应力为对象,其对 应的塑性应变的增量方向即过该应力点的塑性势函 数的法向方向为 r ,主应力空间的塑性修正方向为 u p ? Dr ,对于线性屈服准则而言对应棱线上所有 应力点的塑性修正方向 r p 是相同的, 从而棱线上应 力点的塑性修正方向构成了一个过棱线的边界面, n2?1 的含义见下 如图中向量 l 与 n2?1 组成的面所示, 文定义。

图 6 尖点回映的边界面 Fig.6 Boundary surface of apex

定义尖点处的两个边界面 ③过尖点并且平行于静水压力轴的塑性应变的 增量方向 rt ,对应的塑性修正方向为 u t ? Drt 。 ④过尖点并且平行于压应力平面 ?? 1 ? ? 2 ? 的 塑性应变的增量方向 r1 ,对应的塑性修正方向为
u1 ? Dr1 。

⑤过尖点并且平行于压应力平面 ?? 2 ? ? 3 ? 的 塑性应变的增量方向 r2 ,对应的塑性修正方向为 t u2 ? Dr2 , 则 n1 = u ? u1 , 表示在尖点处由 u t 和 u1 定
区域 2---回映到棱线 区域 1---回映到面
tr 义的边界面的法向量,当 n1 (? ? ? t ) ? 0 时回映到 棱线 l1 ; n2 = u2 ? u t 表示在尖点处由 u t 和 u1 定义的

l
n2 ?1

边界面的法向量,当 n2 (? tr ? ? t ) ? 0 时回映到棱线 l2 。 由于尖点是 l1 和 l2 的交点,可以设定尖点回映 的条件为

up
f =0

n1 ? ? tr ? ? t ? 0 且 n2 ? ? tr ? ? t ? 0

?

?

?

?

(49)

l

5

UMAT 二次开发
ABAQUS 提供了用 FORTRAN 语言编写的子

图 5 边界面示意图 Fig.5 Sketch diagram of boundary surface

程序接口,供用户二次开发之用,UMAT 子程序具 根据试探应力的位置和边界面可以确定应力回 映的区域: (1)落在边界面与屈服面之间的点需要回映到 面上; (2)落在边界面一侧的点需要回映到棱线上。 采用数学方法进行分析:边界面的法向量为 有以下强大的功能
[12-15]



(1)可以定义材料的本构关系,使用 ABAQUS 材料库中没有包含的材料进行计算, 扩充程序功能; (2)几乎可以用于力学行为分析的任何分析过 程,几乎可以把用户材料属性赋予 ABAQUS 中的 任何单元;

第7期

陈培帅等:Hoek-Brown 准则的主应力回映算法及其二次开发

2217

(3)可以和用户子程序“USDFLD”联合使用, 通过“USDFLD”重新定义单元每一物质点上传递到 UMAT 中场变量的数值。 当 ABAQUS 调用用户子程序时,会把当前的 STEP 和 INCREMENT 利用用户子程序的两个实参 KSTEP 和 KINC 传给用户子程序,从而实现了 ABAQUS 主程序与 UMAT 的相互衔接,已知第 n 步结果开始第 n+1 增量步计算的过程如下: (1)初始化 n+1 步的外力向量 Fn ?1 ? Fn ? ?F , 同时初始化节点的位移向量 U n ?1 ? U n ; (2)开始结构的整体平衡迭代,先计算结构的
int 残余应力 Fnres ,其中 F int 为结构内力, ?1 ? Fn ?1 ? F

沿模型垂直方向施加竖向位移,根据模型的轴对称 性质,采取减缩积分的一次轴对称单元(CAX4R) 进行分析,模型如图 8 所示。

然后由一致性刚度矩阵形成整体刚度矩阵 K 求解更 新的节点位移增量 ?U n ?1 ? K ?1F res ,从而计算出应 变的增量 ?? n?1 , 然后开始检查整体平衡迭代的收敛 性:如果 Fnres ?1 ? TOL ,则不收敛,进入第 3 步的计 算;否则进入第 4 步; ( 3) 开始应力更新, 这一步通过 UMAT 子程序 实现,将上一步计算出的应变增量 ?? n?1 传入子程 序,计算出更新后的应力和一致性刚度矩阵,然后 转入第 2 步进行新一轮的迭代; (4)增量步结束,更新 n+1 步的位移 U n ?1 ?
U n ? ?U n ?1 。
图 8 圆柱体有限元模型 Fig.8 Finite element model of cylinder

从图 9~11 的计算结果的可以发现,H-B 模型 与 M-C 模型的计算结果相差不大,不同围压情况 下,H-B 模型的屈服应力与 M-C 模型相比,有 3 种 情况:前者大于后者,如图 9 所示;前者非常接近 于后者,如图 10 所示;前者小于后者,如图 11 所 示。这是因为 H-B 模型为非线性的强度包络线,而 M-C 模型为线性强度包络线,对于同一种岩石材 料,两者的强度包线可能有相交的情况。
100
轴向应力/MPa

6





基于开发的 UMAT 用户子程序, 对室内常规岩 石力学试验进行有限元模拟,并与 ABAQUS 内嵌 的 Mohr-Coulomb 模型的计算结果进行比较,验证 子程序的计算效果。为方便结果对比,本算例的参 数采用非相关联流动法则。如图 7 所示,根据试验 数据,拟合出材料的参数:E = 34.8 GPa,v ? 0.25 ;
a ? 0.5 , s ? 0.151 2 , m ? 2.725 ,? ci ? 230 MPa; c ? 21.569 MPa, ? ? 35.42 °。
300
最大主应力/MPa

80 60 40 20 0 0 1 M-C 模型 H-B 模型 2 3 4 5

-3

-2

-1

轴向应变/‰

侧向应变/‰

图 9 无围压情况下的应力-应变曲线 Fig.9 Stress-strain curves with no confining pressure

250 200
轴向应力/MPa

200 Hoek-Brown 拟合曲线 试验数据 0 10 20 30 40 50 60 -7 -5 150 100 50 0 -3 -1 1 3 5 7 轴向应变/‰ 侧向应变/‰ M-C 模型 H-B 模型

150 100 50 0

-20 -10

最小主应力/MPa

图 7 H-B 模型拟合曲线 Fig.7 Fitting curve of H-B model

基于上述参数,进行了常规三轴压缩试验的数 值模拟,圆柱体高度为 100 mm,直径为 50 mm,

图 10 围压为 30 MPa 情况下的应力-应变曲线 Fig.10 Stress-strain curves with confining pressure of 30 MPa

2218









2011 年

250
轴向应力/MPa

200 150 100 50 0 0 2.5 M-C 模型 H-B 模型

估脆性岩石动态强度的适用性 [J]. 岩石力学与工程学 报, 2003, 22(2): 171-176. ZHAO Jian, LI Hai-bo. Estimating the dynamic strength of rock using Mohr-Coulomb and Hoek-Brown criteria[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(2): 171-176. [6]
10

-7.5 -5.0

-2.5

5.0

7.5

轴向应变/‰

侧向应变/‰

图 11 围压为 40 MPa 情况下的应力-应变曲线 Fig.11 Stress-strain curves with confining pressure of 40 MPa

[7]

7





本文提出了 H-B 准则主应力空间的隐式积分算 法,解决了棱线和尖点处难以收敛的问题,同时也 解决了很多研究中圆滑函数修正法或 D-B 估量法导 致的模型误差的问题。依托 ABAQUS 软件平台, 编制了 UMAT 用户子程序,弥补了 ABAQUS 本构 模型的一个空白,通过常规三轴压缩试验的数值模 拟,对比了 H-B 模型与 M-C 模型的计算结果,表 明 H-B 模型与 M-C 模型在受压情况下的计算结果 比较接近,两者的屈服应力大小有交替现象,验证 了 H-B 模型的非线性性质, 很好地验证了程序的准 确性和有效性。 参 考 文 献
[1] 严春风, 徐健, 陈彦峰, 等. 岩体霍克-布朗经验强度 准则参数随机特征分析[J]. 岩石力学与工程学报, 1999, 18(5): 497-502. YAN Chun-feng, XU Jian, CHEN Yan-feng, et al. Probabilistic analysis on hoek brown criterion parameters for rockmass[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 18(5): 497-502. 吴顺川, 金爱兵, 高永涛. 基于广义霍克-布朗准则的 边坡稳定性强度折减法数值分析 [J]. 岩土工程学报 , 2006, 28(11): 1975-1980. WU Shun-chuan, JIN Ai-bing, GAO Yong-tao. Numerical simulation analysis on strength reduction for slope of jointed rock masses based on generalized Hoek-Brown failure criterion[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(11): 1975-1980. 解联库, 杨小聪, 刘庆林, 等. 基于霍克-布朗破坏准则 的采场空区稳定性分析[J]. 有色金属, 2008, 60(3): 109 -111. XIE Lian-ku, YANG Xiao-cong, LIU Qing-lin. Stability analysis of stope goaf based on Hoek-Brown failure criterion[J]. Nonferrous Metals, 2008, 60(3): 109-111. HOEK E, BROWN E T. Empirical strength criterion for rock masses[J]. Journal of Geotech. Eng., 1980, 106: 1013-1035. 赵坚, 李海波. 摩尔-库仑和霍克-布朗强度准则用于评 [8]

SHARAN S K. Exact and approximate solutions for displacements around circular openings in elastic-brittle-plastic Hoek–Brown rock[J]. Int. J. Rock Mech. Min. Sci., 2005, 42: 542-549. 林杭, 曹平, 李江腾, 等. 基于 Hoek-Brown 准则的三 维边坡变形稳定性分析 [J]. 岩土力学, 2010, 31(11): 3656-3660. LIN Hang, CAO Ping, LI Jiang-teng, et al. Deformation stability of three-dimensional slope based on Hoek-Brown criterion[J]. Rock and Soil Mechanics, 2010, 31(11): 3656-3660.

WAN R G. Implicit integration algorithm for Hoek-Brown elastic-plastic model[J]. Comput. Geotech., 1992, 14: [9] CRISFIELD M A. Non-linear finite element analysis of solids and structures, Vol. 2: advanced topics[M]. New York: John Wiley & Sons, 1997. [10] HOEK E, BROWN E T. Practical estimates of rock mass strength[J]. International Journal of Rock Mechanics and Mining Sciences, 1997, 34(8): 1165-1186. [11] 贾 善 坡 , 陈 卫 忠 , 杨 建 平 , 等 . 基 于 修 正 MohrCoulomb 准则的弹塑性本构模型及其数值实施 [J]. 岩 土力学, 2010, 31(7): 2051-2058. JIA Shan-po, CHEN Wei-zhong, YANG Jian-ping,CHEN Pei-shuai,An elastoplastic constitutive model based on modified Mohr-Coulomb criterion and its numerical implementation[J]. Rock and Soil Mechanics, 2010, 31(7): 2051-2058. [12] JOHAN C,LARS D. An exact implementation of the Hoek–Brown criterion for elasto-plastic finite element calculations[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(6): 831-847. [13] 杨曼娟. ABAQUS 用户材料子程序开发及应用[D]. 武 汉: 华中科技大学, 2005. [14] 范庆来, 栾茂田, 杨庆. 修正剑桥模型的隐式积分算法 在 ABAQUS 中的数值实施[J]. 岩土力学, 2008, 29(1): 269-273. FAN Qing-lai1, LUAN Mao-tian, YANG Qing. Numerical implementation of implicit integration algorithm for modified cam-clay model in ABAQUS[J]. Rock and Soil Mechanics, 2008, 29(1): 269-273. [15] 岑威钧, 朱岳明. 基于 ABAQUS 的土石料本构模型二 次开发及其应用[J]. 水利水电科技进展, 2005, 25(6): 78-81. CEN Wei-jun, ZHU Yue-ming. ABAQUS-based secondary development of constitutive model for earth rockfill materials and its application[J]. Advances in Science and Technology of Water Resources, 2005, 25(6): 78-81.

[2]

[3]

[4]

[5]


相关文章:
更多相关标签: